Could eBird data track population dynamics in time-series in the context of ecological release?

Ecological release: Population level niche expansion and/or shift when antagonistic interspecific interaction is removed or reduced (Herrmann et al., TREE 2021). Hypothetical ecological release could be explored as a function of species richness. Lower competition in islands expand niche (e.g., abundance), while more habitats and more species increase level of specialization in the Mainland (Sherry et al., Auk 2020).

Using the Caribbean meta-archipelago, we filtered the no-aquatic birds reported in eBird, between Julian day 125 and 150 from 2012-2013, between 6 and 11 hours in the morning. We excluded incomplete list and retained checklists of only “traveling” or “stationary” protocol, with a minimum effort of less than 2 km, more than 10 minutes but less than 2 hours (10:120 min) and less than 10 observers. We assumed as counting samples the mean of counts of each species reported by the observers during the 25 days filtered per year; if the mean is very constant among independent observers, we assume a higher probability of detection, while higher variation of the counts might reflect a difference in detection probability. These data constitutes the primary time-series of species for different islands.

To test a model of stochasticity process, we selected a widespread and easy to detect species in the meta-archipelago, which are distributed from the Mainland to different islands in the Caribbean: Bananaquit Coereba flaveola. Theory predicts that species in the mainland have lower density than islands, and that density could increase as species richness and potentially interspecific competition decreases.

Filtering eBird data in meta-archipelagos

The code in HiperGator was (Do not run):

########################################
###~ eBird data - Island filtering ~####
########################################

#Packages
library(auk); #eBird data filters
library(tidyverse); #manipulate data
library(dggridR) #Spatiotemporal Subsampling
library(iNEXT) #rarefaction

#~ eBird data ~####
#Data downloaded from eBird, up to August 2023, saved in the folder of the HPG

#Define the range
CaribbeanRange <- c(-91,8,-59,27)
IndoMalayRange <- c(90,-20,180,20)

#select the columns to extract (based on x$col_idx$name)
colsE <- c("observer_id", "sampling_event_identifier",
           "group identifier",
           "common_name", "scientific_name",
           "observation_count",
           "country", "state_code", "locality_id", "latitude", "longitude",
           "protocol_type", "all_species_reported",
           "observation_date",
           "time_observations_started",
           "duration_minutes", "effort_distance_km",
           "number_observers")

#For Caribbean
f_ebdCa <- "ebd_IslandsCa.txt" #Temporal file to save the filtering eBird data (records)
f_sedCa <- "sed_IslandsCa.txt" #Temporal file to save the filtering sampling event data

ebd_filtCarib<-auk_ebd("eBirddata/ebd_relAug-2023.txt",
                       file_sampling = "eBirddata/ebd_sampling_relAug-2023.txt") %>%
  auk_bbox(CaribbeanRange) %>% #W, S, E, N 
  auk_year(c(2012:2023)) %>%
  auk_protocol(c("Traveling", "Stationary")) %>%
  auk_distance(distance = c(0,5)) %>%
  auk_duration(duration = c(0,300))%>%
  auk_complete() %>%
  auk_filter(f_ebdCa, f_sedCa, overwrite=T, keep = colsE)

f_ebdIM <- "eBirddata/ebd_IslandsIM.txt" #Temporal file to save the filtering eBird data (records)
f_sedIM <- "eBirddata/sed_IslandsIM.txt" #Temporal file to save the filtering sampling event data

ebd_filtIndoP<-auk_ebd("eBirddata/ebd_relAug-2023.txt",
                       file_sampling = "eBirddata/ebd_sampling_relAug-2023.txt") %>%
  auk_bbox(IndoMalayRange) %>% #W, S, E, N 
  auk_year(c(2012:2023)) %>%
  auk_protocol(c("Traveling", "Stationary")) %>%
  auk_distance(distance = c(0,5)) %>%
  auk_duration(duration = c(0,300))%>%
  auk_complete() %>%
  auk_filter(f_ebdIM, f_sedIM, overwrite=T, keep = colsE)

#and with read_ebd I apply another filter to do not repeat records from groups
sed_Carib <- read_sampling(f_sedCa)
ebd_Carib <- read_ebd(f_ebdCa) 
sed_IndoM <- read_sampling(f_sedIM)
ebd_IndoM <- read_ebd(f_ebdIM) 

# Function to convert time observation to hours since midnight
time_to_decimal <- function(x) {
  x <- hms(x, quiet = TRUE)
  hour(x) + minute(x) / 60 + second(x) / 3600
}

#Caribbean data

Carib_eff <- ebd_Carib %>%
  mutate(
    # I don't have here count in X to convert to NA
    observation_count = as.integer(observation_count),
    # effort_distance_km to 0 for non-travelling counts
    effort_distance_km = if_else(protocol_type == "Stationary",
                                 0, effort_distance_km),
    # convert time to decimal hours since midnight
    time_observations_started = time_to_decimal(time_observations_started),
    hour_sampling = round(time_observations_started, 0),
    # split date into year, month, week, and day of year
    year = year(observation_date),
    month = month(observation_date),
    week = week(observation_date),
    day_of_year = yday(observation_date)) %>%
  filter(number_observers <= 10,         #Only list with less than 10 observers
         effort_distance_km <= 5,        #be sure of distance effort
         duration_minutes %in% (0:300), #be sure of duration effort
         !is.na(observation_count))      #only records with abundance estimation

Checklists_with_more_than_fourSppC <- Carib_eff %>%
  select(sampling_event_identifier, scientific_name) %>%
  group_by(sampling_event_identifier) %>%
  summarise(N_species = n()) %>%
  filter(N_species >= 4)
#352797 lists

Caribbean_cooc <- Carib_eff[Carib_eff$sampling_event_identifier %in% Checklists_with_more_than_fourSppC$sampling_event_identifier, ] 

Carib_sed <- sed_Carib %>%
  mutate(
    # effort_distance_km to 0 for non-travelling counts
    effort_distance_km = if_else(protocol_type == "Stationary",
                                 0, effort_distance_km),
    # convert time to decimal hours since midnight
    time_observations_started = time_to_decimal(time_observations_started),
    hour_sampling = round(time_observations_started, 0),
    # split date into year, month, week, and day of year
    year = year(observation_date),
    month = month(observation_date),
    week = week(observation_date),
    day_of_year = yday(observation_date)) %>%
  filter(number_observers <= 10,         #Only list with less than 10 observers
         effort_distance_km <= 5,        #be sure of distance effort
         duration_minutes %in% (0:300)) #be sure of duration effort

#And now with the Oriental-Indo_Malayan-Papuan-Melanesian
IndoM_eff <- ebd_IndoM %>%
  mutate(
    # I don't have here count in X to convert to NA
    observation_count = as.integer(observation_count),
    # effort_distance_km to 0 for non-travelling counts
    effort_distance_km = if_else(protocol_type == "Stationary",
                                 0, effort_distance_km),
    # convert time to decimal hours since midnight
    time_observations_started = time_to_decimal(time_observations_started),
    hour_sampling = round(time_observations_started, 0),
    # split date into year, month, week, and day of year
    year = year(observation_date),
    month = month(observation_date),
    week = week(observation_date),
    day_of_year = yday(observation_date)) %>%
  filter(number_observers <= 10,         #Only list with less than 10 observers
         effort_distance_km <= 5,        #be sure of distance effort
         duration_minutes %in% (0:300), #be sure of duration effort
         !is.na(observation_count))      #remove records without abundance estimation  

Checklists_with_more_than_fourSppIP <- IndoM_eff %>%
  select(sampling_event_identifier, scientific_name) %>%
  group_by(sampling_event_identifier) %>%
  summarise(N_species = n()) %>%
  filter(N_species >= 4)

IndoM_eff_cooc <- IndoM_eff[IndoM_eff$sampling_event_identifier %in% Checklists_with_more_than_fourSppIP$sampling_event_identifier, ]


IndoM_sed <- sed_IndoM %>%
  mutate(
    # effort_distance_km to 0 for non-travelling counts
    effort_distance_km = if_else(protocol_type == "Stationary",
                                 0, effort_distance_km),
    # convert time to decimal hours since midnight
    time_observations_started = time_to_decimal(time_observations_started),
    hour_sampling = round(time_observations_started, 0),
    # split date into year, month, week, and day of year
    year = year(observation_date),
    month = month(observation_date),
    week = week(observation_date),
    day_of_year = yday(observation_date)) %>%
  filter(number_observers <= 10,         #Only list with less than 10 observers
         effort_distance_km <= 5,        #be sure of distance effort
         duration_minutes %in% (0:300)) #be sure of duration effort

#Only those land birds ####
#AVONET filtered by No aquatics
filterAvonet <- read.csv("filterAVONET/AVONET_filter_NoAquatics.csv")

#to combine the data, we have to adjust the names
names(filterAvonet)[2] <- "scientific_name"

#filter in eBird data those species identify in AVONET as land birds
ebird.carib = Caribbean_cooc[Caribbean_cooc$scientific_name %in% filterAvonet$scientific_name, ]
ebird.indop = IndoM_eff_cooc[IndoM_eff_cooc$scientific_name %in% filterAvonet$scientific_name, ]

##Spatial Grids ####
set.seed(123)

# Spatiotemporal Subsampling in diameters of ~1 km (area of ~100 km^2)
dggs_island <- dgconstruct(spacing = 11) #spacing 1 correspond to Characteristic Length Scale, or diameter of spherical cell

#add a new variable that identify cell
Caribbean_SS <- ebird.carib %>%
  mutate(cell = dgGEO_to_SEQNUM(dggs_island, #id for cells
                                longitude, latitude)$seqnum) %>%
  group_by(cell, scientific_name) %>%
  mutate(mu_count = mean(observation_count))

IndoPacific_SS <- ebird.indop %>%
  filter(observation_count < 7000) %>% #there is an outlying count for _Eurylaimus ochromalus_ https://ebird.org/checklist/S112199250
  mutate(cell = dgGEO_to_SEQNUM(dggs_island, #id for cells
                                longitude, latitude)$seqnum) %>%
  group_by(cell, scientific_name) %>%
  mutate(mu_count = mean(observation_count))

#Get the number of species in each cell
RichnessCellCaribbean   <- Caribbean_SS %>% 
  group_by(cell, scientific_name) %>%
  summarise(count=n()) %>% 
  group_by(cell) %>% 
  summarise(SpRichness=n()) %>%
  filter(SpRichness >= 4)

#Exclude cells in SS with less than 4 species
Caribbean_SS = Caribbean_SS[Caribbean_SS$cell %in% RichnessCellCaribbean$cell, ]

#Community assemblages
Caribbean_Cell_Assembly <- Caribbean_SS %>%
  select(cell, scientific_name, mu_count) %>%
  summarise(mu_count = mean(mu_count)) %>%
  pivot_wider(names_from = cell, values_from = mu_count, values_fill = 0)

#Get the grid cell boundaries for cells which had quakes
gridCaribbean <- dgcellstogrid(dggs_island,RichnessCellCaribbean$cell)

#Update the grid cells' properties to include the number of lists in each cell
gridCaribbean <- merge(gridCaribbean,RichnessCellCaribbean,by.x="seqnum",by.y="cell")

# Handle cells that cross 180 degrees
wrapped_gridCaribbean = st_wrap_dateline(gridCaribbean, 
                                         options = c("WRAPDATELINE=YES","DATELINEOFFSET=180"), quiet = TRUE)

#Backup
saveRDS(wrapped_gridCaribbean, "Completeness_data/wrapped_gridCaribbean.rds")
saveRDS(Caribbean_SS, "Completeness_data/Caribbean_SS_140324.rds")
saveRDS(Caribbean_Cell_Assembly, "Completeness_data/Caribbean_Cell_Assembly.rds")

#Get the number of records in each cell
RichnessCellIndoPacific   <- IndoPacific_SS %>% 
  group_by(cell, scientific_name) %>%
  summarise(count=n()) %>% 
  group_by(cell) %>% 
  summarise(SpRichness=n()) %>%
  filter(SpRichness >= 4)

#Exclude cells in SS with less than 4 species
IndoPacific_SS = IndoPacific_SS[IndoPacific_SS$cell %in% RichnessCellIndoPacific$cell, ]

#Community assemblages
IndoPacific_Cell_Assembly <- IndoPacific_SS %>%
  select(cell, scientific_name, mu_count) %>%
  summarise(mu_count = mean(mu_count)) %>%
  pivot_wider(names_from = cell, values_from = mu_count, values_fill = 0)

#Get the grid cell boundaries for cells which had quakes
gridIndoPacific <- dgcellstogrid(dggs_island,RichnessCellIndoPacific$cell)

#Update the grid cells' properties to include the number of lists in each cell
gridIndoPacific <- merge(gridIndoPacific,RichnessCellIndoPacific,by.x="seqnum",by.y="cell")

# Handle cells that cross 180 degrees (and some weird in the map)
wrapped_gridIndoPacific = st_wrap_dateline(gridIndoPacific, 
                                           options = c("WRAPDATELINE=YES","DATELINEOFFSET=180"), quiet = TRUE)

#Backup
saveRDS(wrapped_gridIndoPacific, "Completeness_data/wrapped_gridIndoPacific.rds")
saveRDS(IndoPacific_SS, "Completeness_data/IndoPacific_SS_140324.rds")
saveRDS(IndoPacific_Cell_Assembly, "Completeness_data/IndoPacific_Cell_Assembly.rds")
 
#End of this code

Estimate sample-coverage completeness per sampling unit

Then, we estimated the completeness by sampling units (Cell id in a grid of ~\(100~ km^2\)). This was also conducted in the HiperGator (Do not run):

### Assess completeness of the sampling units

#* Orlando Acevedo-Charry
#* Updated on March 22nd, 2024


#Packages

library(tidyverse)
library(iNEXT) #rarefaction

#Load the data from Code_Islands_eBird_filtering.R

Caribbean_SS <- readRDS("Completeness_data/Caribbean_SS_140324.rds")
IndoPacific_SS <- readRDS("Completeness_data/IndoPacific_SS_140324.rds")

Caribbean_Cell_Assembly <- readRDS("Completeness_data/Caribbean_Cell_Assembly.rds")
IndoPacific_Cell_Assembly <- readRDS("Completeness_data/IndoPacific_Cell_Assembly.rds")

wrapped_gridCaribbean <- readRDS("Completeness_data/wrapped_gridCaribbean.rds")
wrapped_gridIndoPacific <- readRDS("Completeness_data/wrapped_gridIndoPacific.rds")

# Completeness with iNEXT and new grid ####

#Caribbean region
names(Caribbean_Cell_Assembly) #each row is a species (2254), and columns 2:10725 is the cell ID

listCaribbean<-as.data.frame(Caribbean_Cell_Assembly[,2:ncol(Caribbean_Cell_Assembly)])

#iNterpolation and EXTrapolation for the first three Hill numbers (see Chao et al. 2004):
Caribbean_iNEXT<-iNEXT(listCaribbean,q=c(0,1,2), 
                           datatype = "abundance", 
                           endpoint = NULL, #double the reference sample size
                           conf=0.95, nboot=500)
saveRDS(Caribbean_Completeness, "Completeness_output/Caribbean_Completeness.rds") #DONE

CompleteCaribbeanCells <- Caribbean_iNEXT$iNextEst$coverage_based %>%
  filter(Method == "Observed" & Order.q %in% 0 & SC >= 0.9)  

#Exclude cells in SS with less than 4 species
Caribbean_SS = Caribbean_SS[Caribbean_SS$cell %in% CompleteCaribbeanCells$Assemblage, ]

saveRDS(Caribbean_SS, "Completeness_output/Caribbean_SS_Completeness.rds") #DONE

#Get the grid cell boundaries for cells with completeness >0.9
wrapped_gridCaribbean = wrapped_gridCaribbean[wrapped_gridCaribbean$seqnum %in% CompleteCaribbeanCells$Assemblage, ]

gridCaribbean = merge(wrapped_gridCaribbean,CompleteCaribbeanCells,by.x="seqnum",by.y="Assemblage")

saveRDS(gridCaribbean, "Completeness_output/gridCaribbean.rds") #DONE

#Now with the IndoPacific
listIndoPacific<-as.data.frame(IndoPacific_Cell_Assembly[,2:ncol(IndoPacific_Cell_Assembly)])

#iNterpolation and EXTrapolation for the first three Hill numbers (see Chao et al. 2004):
IndoPacific_iNEXT<-iNEXT(listIndoPacific,q=c(0,1,2), 
                         datatype = "abundance", 
                         endpoint = NULL, #double the reference sample size
                         conf=0.95, nboot=500)
saveRDS(IndoPacific_iNEXT, "Completeness_output/IndoPacific_Completeness.rds") #DONE

CompleteIndoPacificCells <- IndoPacific_iNEXT$iNextEst$coverage_based %>%
  filter(Method == "Observed" & Order.q %in% 0 & SC >= 0.9)  

#Exclude cells in SS with less than 4 species
IndoPacific_SS = IndoPacific_SS[IndoPacific_SS$cell %in% CompleteIndoPacificCells$Assemblage, ]

saveRDS(IndoPacific_SS, "Completeness_output/IndoPacific_SS_Completeness.rds") #DONE

#Get the grid cell boundaries for cells with completeness >0.9
wrapped_gridIndoPacific = wrapped_gridIndoPacific[wrapped_gridIndoPacific$seqnum %in% CompleteIndoPacificCells$Assemblage, ]

gridIndoPacific = merge(wrapped_gridIndoPacific,CompleteIndoPacificCells,by.x="seqnum",by.y="Assemblage")
saveRDS(gridIndoPacific, "Completeness_output/gridIndoPacific.rds") #DONE

#End of this code

Example: Bananaquit populations

We selected a widespread distributed species to test our predictions of difference in stationary dynamic abundance from the mainland to islands, the Bananaquit Coereba flaveola

library(tidyverse) #seven packages in one for data management and visualization
library(MASS); #R functions and datasets to support "Modern Applied Statistics with S", a book from W.N. Venables and B.D. Ripley
source("ROUSSE-2.0.R") #Code from Dennis & Ponciano Ecology 2014.

#Call data and manipulate it 

Coefla <- readRDS("Completeness_output/Caribbean_SS_Completeness_elevation.rds") |>
  ungroup() |>
  filter(scientific_name == "Coereba flaveola") |>
  mutate(Day.t = as.numeric(observation_date)-15339, #First estimation of estimate in a year basis
         Site.k = as.character(cell), 
         Year.t = as.numeric(year)) 

#Filter cells with more than 3 years of sampling
CellsMore3yrs <- Coefla %>% 
  group_by(Site.k, Year.t) %>%
  summarise(Lists = n()) %>% 
  group_by(Site.k) %>% 
  summarise(List_Year = n()) %>%
  filter(List_Year >= 3)

#Select only cells with more than 3 years of data
Coefla = Coefla[Coefla$Site.k %in% CellsMore3yrs$Site.k, ]

#Mean count per year per site
Coefla <- Coefla |>
  group_by(Site.k, Year.t) |>
  summarise(Observed.year.y = round(mean(observation_count),0)) |>
  left_join(Coefla, by = c("Site.k","Year.t"))

#Mean count per day per site

Coefla <- Coefla |>
  group_by(Site.k, Day.t) |>
  summarise(Observed.day.y = round(mean(observation_count),0)) |>
  left_join(Coefla, by = c("Site.k","Day.t"))

Coefla |> 
  dplyr::select(state_code, Site.k, Year.t, Day.t, Observed.year.y, Observed.day.y, observation_count) |>
  arrange(state_code) |>
  head()

saveRDS(Coefla, file = "Coereba_flaveola_data.rds")

We can use the map of coast line Mainlan-Island map from Sayer et al. (J. Oper. Oceanogr.)

library(sf); #deal with shape files
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.0     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(maps); #load maps
## 
## Attaching package: 'maps'
## 
## The following object is masked from 'package:purrr':
## 
##     map
library(ggmagnify); # Zoom detail (inset figures)

Coefla <- readRDS("Coereba_flaveola_data.rds")

#The grid of cells
gridCaribbeanPost <- readRDS("Completeness_output/gridCaribbean.rds") #3823 Sites

#Only have the values in Coefla
gridCaribbeanPost = gridCaribbeanPost[gridCaribbeanPost$seqnum %in% Coefla$Site.k, ] #1033 Sites (cells k)

#A world map as reference

world1 <- sf::st_as_sf(map(database = 'world', plot = FALSE, fill = TRUE))
world1
## Simple feature collection with 253 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -85.19218 xmax: 190.2708 ymax: 83.59961
## Geodetic CRS:  +proj=longlat +ellps=clrk66 +no_defs +type=crs
## First 10 features:
##                                        ID                           geom
## Aruba                               Aruba MULTIPOLYGON (((-69.89912 1...
## Afghanistan                   Afghanistan MULTIPOLYGON (((74.89131 37...
## Angola                             Angola MULTIPOLYGON (((23.9665 -10...
## Anguilla                         Anguilla MULTIPOLYGON (((-63.00122 1...
## Albania                           Albania MULTIPOLYGON (((20.06396 42...
## Finland                           Finland MULTIPOLYGON (((20.61133 60...
## Andorra                           Andorra MULTIPOLYGON (((1.706055 42...
## United Arab Emirates United Arab Emirates MULTIPOLYGON (((53.92783 24...
## Argentina                       Argentina MULTIPOLYGON (((-64.54916 -...
## Armenia                           Armenia MULTIPOLYGON (((45.55235 40...
#Global Shoreline Vector
#coast <- st_read("USGSEsriWCMC_GlobalIslands_v3/v10/globalislandsfix.gdb") 

#Different projections!
#Extract CRS of world1 map
world1_crs <- st_crs(world1)
#Transform the geodatabase to match the CRS projection of world1
#coastal <- st_transform(coast, world1_crs)

ggplot() +
  geom_sf(data = world1)+
  geom_sf(data=gridCaribbeanPost,
          aes(color = SpRichness, 
              fill = SpRichness)) +
  scale_color_gradient2(low = alpha("#4575b4", 0.95), 
                        mid = alpha("#ffffbf", 0.95), midpoint = 73,
                        high = alpha("#d73027", 0.95))+
  scale_fill_gradient2(low = alpha("#4575b4", 0.75), 
                       mid = alpha("#ffffbf", 0.75), midpoint = 73,
                       high = alpha("#d73027", 0.75))+
  geom_point(data = Coefla, aes(x = longitude, y = latitude), alpha = 0.5, size = 0.001)+
  coord_sf(xlim = c(-91, -59.4), 
           ylim =  c(8, 27)) +
  scale_x_continuous(breaks = seq(from = -90, to = -60, by = 15))+
  labs(y = "Latitude",
       x = "Longitude",       
       tag = "",
       title = "Neotropical region",
       subtitle = "Caribbean (high eBird coverage-based completeness)",
       color = "Species richness",
       fill = "Species richness") +
  theme_classic()+
  theme(legend.position = "bottom",
        legend.direction = "horizontal")+
  guides(fill = F)+
  geom_magnify(data = Coefla, aes(x = longitude, y = latitude), 
               from = c(xmin = -67.3, xmax = -65.6, ymin = 17.9, ymax = 18.6), 
               to = c(xmin = -70, xmax = -59, ymin = 22, ymax = 26.5),
               shadow = TRUE) 

So, how is the time-series of the observed process?

gridCaribbeanPost$Site.k <- as.character(gridCaribbeanPost$seqnum)

Coefla <- Coefla |>
  left_join(gridCaribbeanPost)
## Joining with `by = join_by(Site.k, seqnum, SpRichness, SC, m, Method, Order.q,
## qD, qD.LCL, qD.UCL, geometry)`
ggplot(Coefla, aes(x = observation_date, group = Site.k))+
  geom_line(aes(y = log(Observed.day.y), color = SpRichness), alpha = 0.25, linetype = "dotted")+
  geom_line(aes(y = log(Observed.year.y), color = SpRichness), alpha = 0.25, linetype = "dashed")+
  scale_color_gradient2(low = alpha("#4575b4", 0.95), 
                        mid = alpha("#ffffbf", 0.95), midpoint = 73,
                        high = alpha("#d73027", 0.95))+
  geom_point(aes(y = log(Observed.day.y)), alpha = 0.01, size = 0.05, shape = 1)+
  geom_point(aes(y = log(Observed.year.y)), alpha = 0.01, size = 0.1, shape = 2)+
  scale_y_continuous(expand = c(0,0))+
  scale_x_date(expand = c(0,0))+
  labs(x = "Observation date",
       y = "Log(mean observed)",
       color = "Species richness")+
  theme_classic()+
  theme(legend.position = "bottom")

This is very messy an difficult to see, lets reduce to values around the Min., 1st Qu., Median, Mean, 3rd Qu., and Max. values of the Species richness gradient.

#Use some examples that represent sp richness:  different species richness

summary(Coefla$SpRichness)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     5.0    37.0    53.0   119.5   225.0   383.0
CoeflaSimp <- Coefla |>
  filter(SpRichness %in% c(12,14,16:29, 31:68,72,
                           104,106,112,115,138,
                           155,156,160,174,179,186:188,196,
                           200:202,205:209,226:231,235,237,243:247,
                           250:253,261,270,278:284,293:299,
                           303:322,330:380)) |> #Select cells with variable richness and more than 200 observations
  group_by(SpRichness) |>
  sample_n(200) |>
  as.data.frame()


ggplot(CoeflaSimp, aes(x = observation_date, group = Site.k))+
  geom_line(aes(y = log(Observed.day.y), color = SpRichness), linetype = "dotted")+
  geom_line(aes(y = log(Observed.year.y), color = SpRichness), linetype = "dashed")+
  scale_color_gradient2(low = alpha("#4575b4", 0.95), 
                        mid = alpha("#ffffbf", 0.95), midpoint = 73,
                        high = alpha("#d73027", 0.95))+
  geom_point(aes(y = log(Observed.day.y)), size = 0.05, shape = 1)+
  geom_point(aes(y = log(Observed.year.y)), size = 0.1, shape = 2)+
  scale_y_continuous(expand = c(0,0))+
  scale_x_date(expand = c(0,0))+
  labs(x = "Observation date",
       y = "Log(mean observed)",
       color = "Species richness")+
  theme_classic()+
  theme(legend.position = "bottom")

In the previous figure, we assumed that the mean \(\mu_{t,k}\) of counts each year (or day in the future) \(t\) in each site \(k\) represents the stationary distribution of the stochastic density dependence in the population.

Look at the lower values of observed abundance in high richness sites, and higher estimates in low richness sites. Could we try to plug-in these time-series in the Gompertz state-state model? The problem is that this data is not uniform… so we need the diffusion process of Ornstein-Uhlenbeck process (Dennis & Ponciano, Ecology 2018).

Lele et al. (Ecol. Lett. 2007) provide the example of the GSS model as a stochastic, density-dependent model (based on Dennis et al., Ecol. Monogr. 2006). In this model, the latent variable of \(N_t\) as the true abundance at time \(t\) in the form: \[ N_t = N_{t-1}* {\rm exp} \left(a + b * {\rm ln}N_{t-1} + E_t \right) \] where \(a\) is the intrinsic growth rate (a constant) and \(b\) is a density dependence parameter (another constant), while the environmental stochasticity is \(E_t \sim {\rm Normal}(0,\sigma^2)\). The variance \(\sigma^2\) is the process noise or environmental variability of the system (Lele et al., Ecol. Lett. 2007), where higher values could serve to identify catastrophic stochasticity (Dennis et al., Ecol. Monogr. 1991).

Thus, we can define \(X_t = {\rm ln}(N_t)\), or the logarithm of the unobserved population abundances at time \(t\), \(t=(0,1,...,q)\). Applying the logarithm scale:

\[ X_t = X_{t-1} + a + (b*X_{t-1}) + E_t \] factorizing out, \[ X_t = a + X_{t-1}(1+b) + E_t \] and simplify the terms, \(c=1+b\) defines the density dependence parameter:

\[ X_t = a + cX_{t-1} + E_t \]

In the GSS model, the probability distribution of \(X_t\) is normal with mean and variance that change in function of time. We begin with a first value \(X_0\), that is plug-in in \(X_1 = a + c(X_0) +E_1\), then this \(X_1\) is plug-in in \(X_2 = a + c(X_1) + E_2 = a + c(a + c(X_0) +E_1) + E_2\), and so on… eventually, when \(X_t \rightarrow{\infty}\), \(X_t\) approaches a time-independent stationary distribution (carrying capacity):

\[ X_{\infty} \sim {\rm Normal} \biggl(\mu=\frac{a}{(1-c)}, {\rm var} = \frac{\sigma^2}{(1-c)^2} \biggr) \] Thus, we can assume that our \(X_0\) is part of that stationary distribution for the long-run time-series.

The observations \(Y_t\) are related to the true population abundances \(X_t\) by: \(Y_t = X_t + F_t\), where \(F_t \sim {\rm Normal}(0,\tau^2)\). Here, the variance \(\tau^2\) quantifies the amount of measurement error or observation error. We have to compute the MLE for \(\Theta = [a, c, \sigma, \tau]\). However, we can also plug the observations \(Y_t\) as a Poisson process (\(Y_t \sim {\rm Poisson}(N_t)\), or in log-scale \(Y_t \sim {\rm Poisson}({\rm ln}(X_t))\)), reducing a parameter estimation (\(\tau\)). However, the nature of the data from eBird require account for the variaton of observation error in the state-space model. In addition, this data is unequal in time intervals, requiring including the diffusion Ornstein-Uhlenbeck process (Dennis & Ponciano, Ecology 2018).

head(CoeflaSimp)
##   Site.k Day.t Observed.day.y Year.t Observed.year.y checklist_id common_name
## 1 743050  1961              3   2017               3    S37017271  Bananaquit
## 2 743050  3667              1   2022               7     G7739403  Bananaquit
## 3 743050  1962              2   2017               3    S37017210  Bananaquit
## 4 598530  1222              2   2015               2    S23327166  Bananaquit
## 5 602177  2913              2   2019               2     G4829124  Bananaquit
## 6 743780  3017              2   2020               5    S66777258  Bananaquit
##    scientific_name observation_count    country state_code locality_id latitude
## 1 Coereba flaveola                 3 Guadeloupe        GP-    L5863751 16.17101
## 2 Coereba flaveola                 1 Guadeloupe        GP-    L2357226 16.16975
## 3 Coereba flaveola                 2 Guadeloupe        GP-    L5863741 16.17204
## 4 Coereba flaveola                 2    Bahamas      BS-EX    L3625605 24.71089
## 5 Coereba flaveola                 2    Bahamas      BS-EX     L831459 24.38150
## 6 Coereba flaveola                 2 Guadeloupe        GP-    L2391989 16.17064
##   longitude observation_date time_observations_started
## 1 -61.12085       2017-05-14                 12.350000
## 2 -61.11943       2022-01-14                 14.216667
## 3 -61.11568       2017-05-15                  8.516667
## 4 -76.82254       2015-05-06                 17.000000
## 5 -76.62037       2019-12-22                 11.000000
## 6 -61.10903       2020-04-04                 12.033333
##                          observer_id        sampling_event_identifier
## 1                         obsr307806                        S37017271
## 2 obsr206146,obsr1979263,obsr2901391 S100736979,S101011323,S101137954
## 3                         obsr307806                        S37017210
## 4                         obsr196942                        S23327166
## 5            obsr1378293,obsr1474099              S62941819,S62990988
## 6                         obsr307806                        S66777258
##   protocol_type duration_minutes effort_distance_km number_observers
## 1    Stationary                5              0.000                1
## 2     Traveling               72              3.800                3
## 3    Stationary               11              0.000                1
## 4     Traveling               30              0.322                1
## 5     Traveling               60              0.600                2
## 6    Stationary                5              0.000                1
##   all_species_reported group_identifier hour_sampling year month week
## 1                 TRUE             <NA>            12 2017     5   20
## 2                 TRUE         G7739403            14 2022     1    2
## 3                 TRUE             <NA>             9 2017     5   20
## 4                 TRUE             <NA>            17 2015     5   18
## 5                 TRUE         G4829124            11 2019    12   51
## 6                 TRUE             <NA>            12 2020     4   14
##   day_of_year   cell mu_count       ID Elevation_STR_90m seqnum SpRichness
## 1         134 743050 3.071429  1328296                 3 743050         12
## 2          14 743050 3.071429  8632992                 3 743050         12
## 3         135 743050 3.071429  1433352                 4 743050         12
## 4         126 598530 3.833333   539992                 6 598530         12
## 5         356 602177 3.457143 10708100                 5 602177         12
## 6          95 743780 5.350000  4546662                 4 743780         12
##          SC  m   Method Order.q qD   qD.LCL   qD.UCL
## 1 1.0000000 34 Observed       0 12 4.174038 19.82596
## 2 1.0000000 34 Observed       0 12 4.174038 19.82596
## 3 1.0000000 34 Observed       0 12 4.174038 19.82596
## 4 0.9227250 25 Observed       0 12 2.481840 21.51816
## 5 1.0000000 33 Observed       0 12 5.755010 18.24499
## 6 0.9266436 26 Observed       0 12 2.673333 21.32667
##                         geometry
## 1 POLYGON ((-61.09389 16.2064...
## 2 POLYGON ((-61.09389 16.2064...
## 3 POLYGON ((-61.09389 16.2064...
## 4 POLYGON ((-76.7697 24.74715...
## 5 POLYGON ((-76.57627 24.4039...
## 6 POLYGON ((-61.01139 16.1688...
CoeflaSimpAnual <- CoeflaSimp |>
  dplyr::select(Site.k, Year.t, Observed.year.y, SpRichness, country) |>
  unique() |>
  filter(Site.k != "2737266")

table(CoeflaSimpAnual$Site.k, CoeflaSimpAnual$Year.t) #this shows that some cells have only 1 or less observations
##          
##           2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023
##   2704578    0    0    0    0    0    1    0    1    0    0    0    0
##   2716904    0    0    0    0    0    0    0    1    0    0    0    0
##   2718358    0    0    0    0    1    1    1    1    1    1    1    1
##   2721275    1    1    1    1    1    1    1    1    1    1    1    1
##   2722003    0    0    0    1    0    0    1    0    1    0    1    1
##   2722005    1    0    0    0    1    1    1    0    1    1    1    1
##   2722726    0    0    0    0    1    0    0    0    0    0    0    0
##   2722732    0    0    1    0    1    1    1    1    1    1    1    1
##   2722733    0    1    1    1    1    1    1    1    1    1    1    1
##   2722734    0    0    1    0    1    1    1    1    1    1    1    1
##   2723462    1    0    0    0    0    1    1    0    0    1    1    1
##   2723463    0    0    0    0    0    0    1    0    0    0    1    1
##   2736545    0    0    0    1    1    0    1    1    0    0    0    0
##   2737996    0    0    0    1    1    2    2    2    2    2    2    1
##   2738030    0    0    0    0    0    0    0    1    0    0    0    0
##   2738759    0    0    0    0    0    0    0    1    0    1    1    0
##   2740215    1    1    1    1    1    1    1    1    1    1    1    1
##   2740216    0    0    0    0    0    0    0    0    0    0    0    1
##   2741648    0    0    0    1    0    1    1    1    1    1    1    0
##   2742392    0    0    0    0    0    1    0    0    0    1    0    0
##   2743848    1    1    0    0    1    1    1    1    1    1    1    1
##   2746771    1    1    1    1    1    1    1    1    0    1    1    1
##   2746772    0    0    0    0    1    1    1    0    1    0    0    0
##   2753407    0    0    1    1    1    1    1    1    1    1    0    0
##   2754136    0    1    1    1    1    1    1    1    1    1    1    1
##   2754137    1    0    0    0    1    0    1    0    0    0    0    1
##   2754865    0    0    0    0    0    0    0    0    0    0    1    1
##   2754866    1    1    1    1    1    1    1    1    1    1    1    1
##   2754867    0    0    1    1    1    1    1    1    1    1    1    1
##   2754868    0    0    1    1    1    0    0    1    1    1    1    1
##   2755597    0    1    0    0    0    0    0    0    1    0    0    1
##   2756337    0    0    0    1    0    1    1    1    1    1    1    1
##   2756338    1    0    0    0    0    1    0    1    0    1    1    0
##   2757068    0    1    0    0    0    1    0    1    1    1    1    1
##   2757069    0    1    1    1    1    1    1    1    1    1    1    1
##   2757070    0    0    0    0    0    1    1    1    1    1    1    1
##   2803707    0    0    0    1    0    1    1    1    1    1    1    1
##   2808766    0    0    0    0    0    0    0    1    1    1    1    0
##   2808806    0    0    1    1    1    1    1    1    1    1    1    1
##   2809535    0    1    0    0    1    1    1    0    1    1    1    1
##   2810226    1    1    1    1    1    1    1    1    1    1    1    1
##   2810230    0    0    1    1    1    1    1    1    1    1    1    1
##   2810956    0    1    1    1    1    1    1    1    1    1    1    1
##   2813875    0    0    1    1    1    1    1    1    1    1    1    1
##   2813881    1    1    1    0    0    1    1    1    1    1    1    1
##   2814603    1    1    1    1    1    1    1    1    1    1    1    1
##   2814605    0    1    1    1    1    1    1    1    1    1    1    1
##   2816069    1    1    1    1    1    1    1    1    1    1    1    1
##   2816798    1    1    1    1    1    1    1    1    1    1    1    1
##   2816799    0    0    1    1    1    1    1    1    1    1    1    1
##   2817526    1    1    1    1    1    1    1    1    1    1    1    1
##   2818258    0    1    0    1    1    1    1    1    1    1    1    1
##   2818987    0    0    0    1    1    1    1    1    1    1    1    1
##   2819718    0    0    0    1    1    1    1    1    1    1    1    1
##   2820437    1    0    0    0    0    0    0    0    0    0    0    1
##   2820448    0    0    0    0    0    1    1    1    1    1    1    1
##   2821166    0    0    0    1    0    1    0    1    0    1    1    1
##   2821168    0    0    0    0    1    1    1    1    1    1    1    1
##   2822633    0    0    0    0    0    1    1    1    1    1    1    1
##   2823363    0    1    0    1    1    1    1    1    1    1    1    1
##   2824088    1    1    1    1    1    1    1    1    1    1    1    1
##   2824094    1    1    1    1    1    1    1    1    1    1    1    1
##   2825550    0    0    0    0    0    0    0    1    1    0    1    1
##   2826280    0    0    1    1    1    1    1    1    1    1    1    1
##   2827736    1    1    1    1    1    1    1    1    1    1    1    1
##   2827747    0    1    1    1    1    1    1    1    1    1    1    1
##   2828477    0    0    1    1    0    1    1    1    1    1    1    1
##   2829207    1    0    1    1    1    1    1    1    1    1    1    1
##   2829937    0    0    0    0    0    0    0    0    2    2    2    2
##   2832840    0    0    0    0    1    1    1    1    1    1    1    1
##   2832841    0    1    0    1    1    1    1    1    1    1    1    1
##   2833572    1    1    0    0    1    1    1    1    1    1    1    1
##   2833573    0    0    0    1    1    1    1    1    1    1    1    1
##   2834299    0    0    0    1    1    1    1    1    1    1    1    1
##   2834300    0    0    0    1    1    1    1    1    1    1    1    1
##   2834304    1    1    0    1    1    1    1    1    1    1    1    1
##   2834306    1    1    1    1    1    1    1    1    1    1    1    1
##   2834307    1    1    0    1    1    1    1    1    1    1    1    1
##   2835030    1    1    1    1    1    1    1    1    1    1    1    1
##   2835031    1    0    1    1    1    1    1    1    1    1    1    1
##   2835036    1    0    0    1    1    1    1    1    1    1    1    1
##   2835764    0    0    1    1    1    1    1    1    1    1    1    1
##   2836492    0    1    0    0    0    0    0    1    1    1    1    1
##   2836493    0    0    1    1    1    1    1    1    1    1    1    1
##   2837953    0    0    0    0    0    0    1    1    1    1    1    1
##   2840143    0    0    0    1    0    1    1    1    1    1    1    1
##   2843066    0    0    0    0    0    0    0    0    0    0    0    1
##   2843088    1    1    1    1    1    1    1    1    1    1    1    1
##   2845278    0    0    0    1    0    0    0    0    1    0    0    1
##   2846008    1    0    1    1    1    1    0    1    0    1    1    1
##   2847460    1    1    1    1    1    1    1    1    1    1    1    1
##   2847461    1    1    1    1    1    1    1    1    1    1    1    1
##   576642     0    1    0    1    1    0    0    0    0    0    0    0
##   577371     0    0    0    0    0    0    1    0    0    0    0    0
##   577372     0    0    1    0    1    1    1    1    0    0    1    0
##   577374     0    0    0    1    0    0    0    1    0    0    0    0
##   579571     0    0    1    1    0    1    1    0    0    0    0    1
##   580300     0    0    1    1    1    1    1    1    0    0    0    0
##   581030     0    0    0    0    1    0    1    0    0    0    0    0
##   581761     0    0    0    1    1    1    1    1    0    0    0    0
##   582490     0    1    0    0    1    1    1    1    0    0    0    1
##   583218     1    0    0    1    1    1    1    1    0    0    0    0
##   583946     0    0    0    0    0    0    0    0    1    0    0    0
##   583947     0    0    0    0    0    0    0    0    0    1    1    0
##   583948     0    0    0    0    0    1    0    0    1    0    1    0
##   584672     0    0    0    0    0    1    0    1    0    0    1    1
##   584675     0    0    0    0    0    0    0    1    0    0    0    0
##   584676     0    0    1    0    1    0    1    1    0    1    0    1
##   585397     0    0    0    0    0    1    1    1    0    1    1    1
##   585402     1    0    0    1    1    0    1    1    0    0    1    1
##   585404     0    0    1    0    0    0    0    0    1    0    1    0
##   586133     0    0    0    0    1    0    0    0    0    0    0    0
##   590499     0    0    0    0    0    1    0    0    0    0    0    0
##   591228     0    1    0    0    1    1    0    0    0    0    0    0
##   591970     0    0    0    0    1    0    0    0    0    0    0    0
##   592686     0    1    0    0    0    0    1    0    0    0    0    0
##   592700     0    0    0    0    0    0    0    0    0    1    1    1
##   593416     0    0    0    0    0    0    0    0    0    1    0    1
##   593420     0    0    0    0    0    0    0    0    1    0    0    0
##   593421     0    0    0    0    0    0    0    1    0    0    0    0
##   594149     0    0    1    0    1    0    0    1    1    1    0    0
##   594150     0    0    1    0    1    0    0    0    0    1    1    1
##   594151     0    1    1    0    1    1    1    1    1    1    1    1
##   594152     0    0    0    0    0    0    0    1    0    1    0    0
##   594159     0    0    0    1    0    0    0    1    0    0    0    1
##   594875     0    0    0    0    0    0    0    0    0    0    1    1
##   594880     0    0    0    0    0    0    0    0    1    0    0    0
##   594881     0    0    0    0    0    0    0    1    0    0    0    0
##   594889     0    1    0    0    0    0    1    0    0    1    1    1
##   595604     0    0    0    0    0    0    0    0    0    0    1    0
##   595605     0    0    0    0    0    0    0    0    1    0    0    0
##   596349     0    0    0    0    0    1    0    0    0    0    0    0
##   597078     0    0    0    1    0    1    0    1    1    0    0    0
##   597079     0    1    0    0    0    0    1    1    1    1    1    0
##   598530     0    0    0    1    0    0    1    0    0    0    0    1
##   598537     0    0    0    1    0    1    0    0    0    0    1    1
##   599265     0    0    0    0    1    0    1    1    1    1    1    1
##   599266     0    0    0    1    0    0    0    0    1    0    0    0
##   599980     0    0    0    0    0    0    1    0    0    0    1    0
##   600709     0    0    0    0    1    0    0    0    0    0    0    0
##   601439     0    0    0    0    0    0    0    0    0    1    1    1
##   601452     0    1    0    1    0    0    1    1    0    0    1    0
##   602177     1    0    0    1    1    0    1    1    1    1    1    1
##   602182     0    0    0    0    0    1    0    0    0    0    0    1
##   602913     0    1    1    1    0    1    1    1    0    0    1    1
##   602915     0    0    0    0    0    0    0    1    0    0    0    1
##   605095     1    0    0    0    0    0    0    0    0    1    1    1
##   610932     0    0    0    0    0    0    1    1    0    0    1    1
##   611675     0    1    0    1    1    1    1    1    1    0    1    1
##   611676     0    1    0    1    1    0    1    1    1    1    1    1
##   612404     0    1    0    0    0    0    1    1    1    0    1    0
##   613121     0    0    0    0    0    0    0    1    1    0    0    1
##   613133     1    0    0    1    0    0    1    1    1    1    1    1
##   613134     0    1    0    0    0    1    1    1    1    0    1    0
##   618961     0    0    1    1    0    0    1    1    0    0    1    0
##   638671     0    0    0    1    0    0    0    1    0    1    1    1
##   638674     0    0    0    1    0    0    0    0    0    0    1    0
##   639401     1    1    0    0    0    1    0    1    1    1    1    1
##   639404     0    1    1    1    1    1    0    1    0    1    1    1
##   640130     0    1    0    0    0    1    0    0    0    1    1    1
##   640131     0    0    1    1    1    1    1    1    0    1    1    1
##   640863     0    0    1    1    0    0    1    1    1    0    1    1
##   641574     0    0    0    1    0    1    0    0    0    0    0    0
##   641575     0    0    0    0    1    1    1    0    0    0    0    0
##   642303     0    0    0    1    1    1    1    1    1    1    1    1
##   642306     0    0    0    1    0    1    0    0    0    0    0    0
##   643033     0    0    0    1    1    1    0    1    0    1    1    0
##   647358     0    0    0    0    0    1    1    0    0    1    1    1
##   648087     0    1    1    1    0    1    1    1    1    0    1    1
##   648088     1    0    0    0    0    0    0    0    1    0    1    0
##   648090     0    0    0    0    0    0    1    0    0    1    0    0
##   648091     0    0    0    0    0    0    0    1    1    0    1    0
##   648092     1    0    0    0    1    1    1    1    1    1    1    1
##   648817     0    0    0    0    0    0    1    0    0    0    0    0
##   648820     0    1    0    0    0    0    0    0    0    0    0    0
##   648821     1    1    0    1    1    1    1    1    0    1    1    1
##   648823     0    0    0    0    0    0    0    0    0    0    1    0
##   649553     0    0    1    0    1    0    0    0    1    0    1    1
##   649554     0    0    0    0    0    1    1    1    1    0    0    0
##   650278     0    0    0    0    0    0    0    1    0    0    0    0
##   650285     0    0    1    1    0    0    1    1    0    1    1    1
##   651007     0    1    1    0    1    0    1    0    1    1    1    1
##   651010     0    0    0    1    0    1    0    0    0    0    0    0
##   651013     0    0    0    1    1    1    1    1    1    0    0    1
##   651014     0    0    0    0    0    0    0    0    1    1    1    0
##   651015     0    0    0    0    1    0    0    0    0    0    1    0
##   651737     0    1    0    0    0    1    0    0    1    0    0    0
##   651739     0    0    0    0    1    0    0    0    0    0    0    0
##   651745     0    0    0    0    0    1    1    1    0    0    0    0
##   651746     0    1    1    0    0    0    1    0    0    1    1    0
##   652467     0    0    0    0    0    1    0    0    1    0    0    0
##   652468     0    0    1    0    0    0    0    1    1    0    0    0
##   652475     0    0    0    0    0    0    0    1    1    0    0    1
##   653199     0    1    0    1    1    0    0    0    1    0    1    1
##   653936     0    0    0    1    1    1    0    1    1    1    1    0
##   654665     0    0    0    0    0    0    0    0    1    0    0    0
##   655394     0    0    0    0    0    1    0    1    0    0    1    1
##   655395     0    0    0    0    0    0    0    1    0    0    0    0
##   656120     0    0    0    0    0    0    0    0    1    0    0    1
##   656122     0    0    0    0    0    0    0    0    0    0    1    0
##   656123     1    1    0    1    1    1    1    1    1    1    1    1
##   656124     0    0    1    1    1    1    1    1    1    1    1    1
##   656126     0    0    0    1    0    0    0    0    0    1    1    1
##   656850     0    0    0    0    0    0    0    1    0    0    0    0
##   656851     0    0    0    0    0    0    1    0    0    0    1    0
##   656852     0    0    0    1    0    1    1    0    1    1    1    1
##   656853     0    0    0    0    0    0    0    1    1    1    1    1
##   656854     0    0    0    0    0    0    0    0    1    0    0    0
##   656856     0    1    1    0    0    1    0    1    0    0    1    1
##   657577     0    0    1    0    0    1    0    1    0    0    0    0
##   657580     0    0    0    0    0    0    0    0    0    0    0    1
##   657585     0    0    0    1    0    1    1    1    1    0    0    1
##   657586     0    0    0    0    0    0    1    0    0    0    1    1
##   657626     0    0    0    0    0    0    1    1    1    0    1    0
##   657627     0    0    0    0    0    1    0    1    0    0    1    1
##   657628     0    0    0    0    0    0    0    1    0    0    0    0
##   658355     0    0    0    0    0    0    1    0    0    0    0    0
##   658356     0    0    0    0    0    0    1    1    1    0    1    0
##   659043     0    0    0    0    0    0    0    0    0    0    0    1
##   659086     0    0    0    0    0    0    0    1    0    0    0    0
##   659091     0    0    0    0    0    0    1    0    0    1    0    1
##   659820     0    0    0    1    0    0    0    0    1    0    0    0
##   659821     0    0    0    0    0    1    0    0    0    0    0    0
##   660549     0    0    1    1    0    0    1    0    0    0    0    1
##   660554     0    0    1    0    0    0    0    0    0    0    0    0
##   661285     0    0    0    0    0    1    0    0    1    0    0    0
##   662017     0    0    0    0    0    0    1    0    0    0    1    0
##   662746     1    0    0    1    0    1    0    0    0    1    1    1
##   663455     0    0    0    0    0    0    0    0    0    1    0    0
##   663476     0    1    0    0    0    0    1    1    1    0    1    1
##   663477     0    0    1    0    1    1    1    1    1    1    0    1
##   663478     0    0    1    0    0    0    0    1    1    1    1    0
##   664207     0    0    0    0    0    0    0    0    0    1    1    0
##   664933     0    0    0    0    0    0    0    0    0    1    0    0
##   664936     0    0    0    0    0    0    1    0    0    1    0    1
##   665632     0    0    0    0    0    0    1    0    0    0    0    0
##   665662     1    1    0    0    1    0    0    0    0    1    0    1
##   665663     0    0    0    0    0    0    0    1    0    0    0    1
##   665668     0    0    0    0    0    0    0    0    0    1    0    0
##   666361     1    1    1    1    1    1    1    1    1    0    1    0
##   666392     0    0    0    0    0    0    0    0    1    1    0    1
##   667093     0    0    1    1    1    0    0    1    0    0    1    0
##   667123     0    0    0    0    0    0    0    0    0    1    0    1
##   667125     0    1    0    0    0    0    0    0    0    0    0    1
##   667128     0    0    0    0    0    0    0    0    0    1    1    0
##   667854     0    0    0    0    0    0    0    1    1    1    1    1
##   668563     1    0    1    1    1    1    1    1    1    1    1    0
##   668577     0    0    0    0    0    0    0    0    0    0    0    1
##   668579     0    0    0    0    0    0    0    0    0    0    0    1
##   668584     0    0    0    0    1    0    0    0    0    1    0    0
##   668587     0    0    0    0    0    0    0    0    0    0    1    0
##   669289     0    0    0    0    1    1    1    1    0    0    0    0
##   669292     0    1    0    1    0    1    1    1    1    0    0    0
##   669297     0    0    0    0    0    0    0    0    0    1    0    0
##   669308     0    0    0    0    0    1    1    1    0    0    0    0
##   669309     0    0    0    1    1    1    0    1    1    0    1    1
##   670021     1    0    0    0    1    0    1    1    1    1    0    0
##   670037     0    0    0    0    0    0    0    1    1    0    0    0
##   670038     0    0    0    0    0    0    0    0    0    1    0    0
##   670039     0    0    0    0    0    0    0    0    0    0    1    0
##   670040     0    0    0    0    0    0    0    0    0    0    0    1
##   670048     0    0    0    0    0    0    0    0    1    0    0    1
##   670748     0    0    0    0    1    1    1    0    0    0    0    0
##   670749     1    0    0    0    0    0    1    0    0    0    0    0
##   670750     0    0    1    0    1    1    1    1    0    1    1    0
##   670766     0    0    0    0    0    0    0    0    0    0    1    0
##   670776     0    0    0    0    0    0    1    0    0    0    0    0
##   670778     1    0    0    0    1    0    0    1    0    1    0    0
##   671495     0    0    0    1    0    0    0    0    1    1    0    0
##   671497     0    0    0    0    0    0    0    0    1    0    1    0
##   672214     0    0    0    0    0    0    0    0    0    1    0    0
##   672215     0    0    0    0    0    0    0    0    0    0    0    1
##   672227     0    0    0    0    0    0    1    0    1    1    0    0
##   672237     0    0    0    0    0    0    1    0    0    0    0    0
##   672238     0    0    0    0    0    0    0    0    1    0    0    0
##   672942     0    0    0    0    0    0    0    1    0    1    1    0
##   672943     1    1    1    0    0    1    0    0    0    1    1    1
##   672944     1    1    1    1    1    1    1    1    1    1    1    1
##   672945     1    0    0    0    0    0    0    0    0    0    1    0
##   672965     0    1    1    0    1    1    1    1    1    1    1    1
##   672966     0    0    0    0    0    0    0    0    0    1    0    1
##   673672     0    0    0    0    0    0    0    0    0    0    0    1
##   674401     0    0    0    0    0    0    0    1    0    0    0    0
##   674402     0    0    0    1    0    0    0    0    0    0    1    0
##   674405     0    1    0    1    0    0    0    0    0    0    1    0
##   674406     0    1    0    0    0    0    0    0    0    0    0    0
##   674420     0    0    1    0    0    0    0    0    0    0    1    1
##   675130     0    1    0    0    0    1    1    0    0    0    1    1
##   675131     0    0    0    0    0    0    0    0    1    0    0    0
##   675864     0    0    0    0    0    0    0    0    1    1    1    1
##   675877     0    0    0    0    0    0    0    1    1    1    0    0
##   676589     0    0    0    0    0    1    0    0    0    0    0    0
##   676593     0    0    0    0    0    0    0    0    0    0    1    0
##   676594     0    1    0    0    0    0    0    0    1    0    1    0
##   676605     0    0    0    0    0    0    0    0    0    0    0    1
##   676606     1    1    1    1    1    0    0    0    1    1    1    1
##   676608     0    0    0    0    0    1    0    0    0    0    0    0
##   677323     0    0    0    0    0    0    0    1    1    0    0    0
##   677329     0    0    0    0    0    0    0    0    1    1    1    1
##   677334     0    0    0    0    0    0    0    0    1    1    0    0
##   677335     0    0    0    0    0    0    1    1    1    1    1    1
##   677336     1    0    1    1    0    1    1    1    1    0    1    1
##   678049     0    0    0    0    0    0    0    0    1    0    0    0
##   678779     0    0    0    0    0    1    0    0    0    0    0    0
##   678801     0    0    0    1    0    0    0    0    0    0    0    0
##   679530     0    0    1    0    0    0    0    0    0    0    0    0
##   680261     0    0    0    1    0    0    0    0    1    0    0    1
##   680266     0    0    0    0    0    0    0    0    0    0    1    1
##   680991     0    0    0    1    1    0    0    0    1    0    1    1
##   680996     0    0    0    0    0    0    0    0    0    0    0    1
##   681720     0    1    0    0    0    1    1    1    1    1    1    1
##   681725     0    1    1    1    1    1    1    1    1    1    1    1
##   682455     0    1    0    1    1    1    1    1    1    1    1    1
##   683180     0    0    0    0    0    0    1    1    0    1    0    1
##   683181     1    0    0    0    0    0    0    0    0    0    0    0
##   683184     1    1    1    1    1    1    1    1    1    1    1    1
##   683909     0    0    0    0    0    1    0    0    0    0    0    0
##   691218     0    0    1    0    0    1    1    1    1    1    1    1
##   691946     0    1    1    1    1    1    1    1    1    1    1    1
##   691947     0    0    0    0    0    0    0    0    0    0    1    0
##   691948     1    0    1    0    1    1    1    1    1    1    1    1
##   691949     0    0    0    1    1    0    0    1    1    1    1    0
##   692676     1    1    1    1    1    1    1    1    1    0    1    0
##   692678     1    1    1    0    1    1    1    1    1    0    0    0
##   692679     0    0    0    0    1    1    1    1    1    1    1    1
##   693405     0    0    1    1    0    1    1    1    1    1    1    1
##   693406     0    1    0    0    0    0    1    1    1    1    1    0
##   693407     1    0    0    1    1    0    1    1    1    1    1    0
##   693408     0    1    0    0    0    1    1    1    1    1    1    1
##   693409     0    1    1    1    1    1    1    1    1    1    1    1
##   694133     0    1    0    0    1    0    1    0    1    1    1    0
##   694134     1    0    0    1    1    0    0    1    1    0    0    0
##   694135     0    0    0    0    0    0    0    1    0    1    1    0
##   694136     0    0    0    0    1    0    0    1    1    1    1    1
##   694137     0    0    0    0    0    0    0    1    1    0    0    1
##   694138     0    0    0    0    0    0    0    1    0    1    1    0
##   694139     1    0    0    0    0    1    1    1    1    1    1    1
##   694140     1    0    0    0    1    0    0    1    1    1    1    1
##   694862     0    0    1    1    1    0    0    1    1    1    1    1
##   694863     1    1    1    1    1    1    1    1    1    1    1    1
##   694864     0    0    0    0    0    0    0    0    0    1    0    1
##   694865     1    1    1    1    1    1    1    1    1    1    1    1
##   694866     0    0    0    1    0    0    1    1    0    1    1    1
##   694867     0    1    0    1    1    0    0    0    0    1    0    1
##   694868     0    0    0    0    0    0    0    0    0    1    0    0
##   694869     0    1    1    1    1    1    1    1    1    0    1    0
##   694870     0    1    1    1    1    1    1    0    1    1    1    1
##   695592     0    0    1    0    0    0    0    0    1    1    0    1
##   695593     1    1    1    1    1    1    1    1    1    1    1    1
##   695594     0    0    0    1    1    0    0    1    1    0    1    1
##   695595     0    1    0    0    1    0    1    0    0    1    0    0
##   695596     0    0    1    1    0    1    0    0    0    1    0    0
##   695597     0    0    0    1    0    0    0    0    0    0    0    0
##   695598     1    0    1    1    1    1    1    1    1    1    1    1
##   695599     0    0    0    0    1    0    1    1    1    1    1    0
##   695600     1    1    1    1    1    1    1    1    1    1    1    1
##   695601     0    1    1    0    0    0    0    1    1    0    1    1
##   696323     0    1    0    0    1    0    1    0    0    0    1    1
##   696325     1    0    0    0    1    0    0    0    1    0    0    0
##   696326     0    0    0    0    1    0    1    1    1    1    1    1
##   696327     0    0    0    0    0    0    0    1    1    0    1    1
##   696328     1    0    1    1    1    1    1    0    0    1    1    1
##   696329     1    1    1    1    0    0    1    1    1    1    1    0
##   696330     1    1    1    1    0    1    1    1    1    1    1    1
##   696331     0    1    0    0    0    1    0    1    0    1    1    1
##   697053     0    0    0    1    0    0    0    1    0    1    1    0
##   697054     1    1    1    1    1    1    1    1    1    1    1    1
##   697055     1    1    1    1    1    1    0    0    1    1    0    0
##   697056     1    1    1    1    1    0    1    0    0    1    1    1
##   697057     0    0    0    0    0    0    0    1    0    0    1    1
##   697058     0    0    0    1    0    0    0    1    0    0    0    0
##   697059     0    1    1    0    0    1    0    1    0    1    1    1
##   697060     0    0    0    0    0    0    0    0    1    1    1    0
##   697061     1    0    1    1    1    1    1    1    1    1    1    1
##   697784     0    1    1    0    0    0    0    0    0    0    0    0
##   697785     0    1    1    0    0    1    1    1    1    1    1    1
##   697786     0    1    1    0    1    1    0    0    1    1    1    1
##   697787     0    0    0    0    0    0    0    0    1    0    0    1
##   697788     0    0    0    1    0    0    0    1    1    0    1    0
##   697789     1    0    0    0    0    0    0    0    0    0    1    0
##   697790     0    0    0    0    1    0    0    0    0    0    1    0
##   697791     0    0    0    0    1    0    1    1    1    1    1    1
##   697792     1    1    1    1    1    1    1    1    1    1    1    1
##   698515     0    1    0    0    1    1    0    1    1    0    1    1
##   698516     0    0    0    0    1    1    1    0    1    1    0    1
##   698517     0    0    0    0    0    0    1    0    0    0    0    1
##   698518     1    0    0    1    1    0    1    1    0    1    1    0
##   698519     0    0    0    1    1    0    0    1    0    1    1    1
##   698520     0    0    0    0    0    1    0    0    1    1    1    1
##   698521     1    1    1    1    1    1    1    1    1    1    1    1
##   698522     0    1    0    0    0    1    1    1    1    1    1    1
##   699245     0    0    0    0    0    0    0    0    1    0    0    1
##   699247     0    0    1    0    0    0    0    1    1    1    1    0
##   699248     0    0    0    1    1    0    0    0    1    1    1    1
##   699249     1    0    1    1    1    0    1    1    0    1    1    1
##   699250     0    0    0    0    0    0    0    1    1    1    1    1
##   699251     1    1    1    1    0    1    1    1    1    1    1    1
##   699252     1    1    1    1    1    1    1    1    1    1    1    1
##   699976     1    0    0    0    0    0    1    0    0    0    0    0
##   699977     1    0    0    0    0    0    1    0    0    0    1    0
##   699978     0    0    0    0    0    1    0    0    1    1    1    1
##   699979     0    0    0    1    0    0    0    0    0    1    0    1
##   699980     0    1    0    0    1    1    1    1    1    1    1    1
##   699981     0    0    0    0    0    0    1    0    0    0    1    1
##   699982     1    1    0    0    1    0    1    1    1    1    1    1
##   699983     0    0    0    0    0    0    0    0    0    0    0    1
##   700706     0    0    0    0    1    0    1    1    1    1    1    1
##   700707     0    0    0    0    0    0    1    0    0    1    1    1
##   700708     0    1    1    0    0    1    1    1    1    1    1    1
##   700709     0    0    0    0    0    0    1    0    1    1    1    1
##   700710     0    0    0    1    0    0    0    1    0    1    1    1
##   700711     0    0    0    1    1    1    1    0    1    1    1    1
##   700712     1    1    1    1    1    1    1    1    1    1    1    1
##   700713     0    0    0    0    0    0    0    0    0    1    1    1
##   701436     1    1    1    1    1    1    1    1    1    1    1    1
##   701437     0    0    0    0    0    1    1    1    0    1    1    1
##   701438     0    0    0    0    1    0    0    0    1    1    1    1
##   701439     1    0    0    0    1    0    1    0    0    0    1    0
##   701440     0    0    0    0    0    0    0    1    1    0    1    0
##   701441     0    1    1    1    1    1    1    1    1    1    1    1
##   701442     0    0    0    1    0    1    1    1    0    1    1    1
##   701443     0    0    0    0    1    0    0    1    1    0    1    1
##   702167     0    0    0    0    0    0    0    0    0    1    0    0
##   702169     0    0    0    0    0    0    0    0    0    1    1    0
##   702170     1    1    1    1    1    1    1    1    1    1    1    1
##   702171     0    0    0    0    0    0    0    0    0    1    1    0
##   702172     0    1    1    0    1    1    0    1    1    1    1    1
##   702898     0    1    1    1    0    0    0    0    1    1    1    0
##   702899     0    0    0    1    0    0    0    0    0    0    0    1
##   702901     0    0    0    0    0    1    0    1    1    1    1    1
##   703630     1    0    0    0    0    1    0    0    1    1    1    1
##   703631     0    0    0    1    0    1    0    0    0    1    1    1
##   703633     0    1    0    0    0    0    1    0    0    0    0    0
##   703634     0    0    0    1    1    1    1    1    0    0    1    1
##   704360     0    0    0    0    0    0    0    0    0    1    1    1
##   704361     1    1    1    0    1    1    0    1    1    1    1    1
##   704363     0    0    1    1    0    1    1    1    0    0    1    1
##   705091     0    0    0    1    0    0    0    1    0    1    1    0
##   705095     0    0    0    0    0    0    0    1    0    0    0    0
##   705825     1    1    1    1    1    1    1    1    1    1    1    1
##   705826     0    0    0    0    0    0    0    1    1    1    1    1
##   705827     0    0    0    0    1    0    0    1    0    0    1    1
##   705832     0    0    0    1    0    1    0    0    0    0    1    1
##   706555     0    1    1    1    1    0    0    1    0    1    1    1
##   706556     1    1    1    1    1    1    1    1    1    1    1    1
##   706557     1    0    1    0    0    0    0    0    0    0    0    0
##   707286     1    1    1    1    1    1    1    1    1    1    1    1
##   707287     1    0    1    0    1    0    0    0    0    0    1    1
##   707288     1    1    0    0    1    1    0    0    0    0    1    1
##   708017     0    0    0    0    0    0    0    0    0    1    0    1
##   708018     0    1    0    1    1    0    0    1    1    0    1    1
##   708019     1    0    0    0    0    0    0    1    0    0    0    0
##   710200     0    0    0    0    0    0    1    1    0    0    0    1
##   710929     0    0    0    1    1    1    1    1    1    1    1    1
##   710930     0    0    0    1    1    1    1    1    1    1    1    1
##   711659     0    1    1    1    1    1    0    0    1    0    0    1
##   711660     1    1    1    0    1    1    1    1    1    1    1    1
##   711661     1    0    0    1    1    1    1    1    1    1    1    1
##   712390     0    0    1    0    0    1    1    1    1    1    1    1
##   717511     1    1    1    1    1    1    1    1    1    1    1    1
##   717512     0    1    1    1    1    1    1    1    1    1    1    1
##   718240     0    0    1    0    0    1    1    2    0    2    1    2
##   718241     0    1    1    1    1    1    0    1    1    1    1    1
##   718970     1    1    2    1    1    1    1    2    1    1    2    2
##   720425     0    0    0    0    1    0    0    1    0    0    0    0
##   720430     0    0    0    0    0    1    0    1    1    1    1    1
##   721155     0    1    0    0    1    0    0    0    1    0    0    0
##   721160     0    0    0    0    0    0    0    0    1    1    1    1
##   723345     0    0    1    0    1    1    1    0    1    1    1    0
##   723942     0    0    0    0    1    1    1    1    1    0    1    1
##   724679     0    0    1    1    0    1    1    1    1    1    1    1
##   724680     0    0    0    0    1    1    1    1    0    1    1    1
##   724703     0    0    0    0    0    0    0    0    1    0    0    0
##   725360     0    0    1    1    1    1    1    1    1    1    1    1
##   725534     0    0    0    0    0    1    1    1    0    1    1    1
##   726124     0    0    0    0    0    1    1    0    0    1    1    1
##   726137     0    0    0    0    1    0    0    0    0    0    0    1
##   726264     0    1    1    1    1    0    1    1    1    1    1    1
##   726853     0    0    0    0    0    0    0    1    0    1    1    1
##   726873     0    0    0    0    0    1    1    1    1    1    1    1
##   726878     1    1    1    1    1    1    1    1    1    1    1    1
##   726885     0    0    0    0    0    0    0    1    0    0    0    0
##   726994     0    1    1    0    1    0    0    1    0    1    1    1
##   727604     0    0    0    0    0    1    0    1    0    0    1    0
##   727605     0    1    0    0    0    0    0    0    0    0    0    1
##   727609     0    0    0    0    0    0    0    0    0    0    0    1
##   727626     0    0    0    0    0    0    1    1    1    0    0    0
##   727723     0    1    0    0    1    0    0    1    1    1    0    1
##   728453     0    0    0    0    1    1    1    1    0    1    0    0
##   728461     0    0    0    0    0    0    0    0    0    1    0    0
##   728462     0    1    0    0    0    1    0    0    1    0    1    1
##   729098     1    1    1    1    1    1    1    1    1    1    1    1
##   729191     0    0    0    0    0    0    1    0    1    0    1    1
##   729827     0    0    1    1    1    1    1    1    0    1    1    1
##   729828     0    1    0    0    0    1    1    1    1    1    1    1
##   730557     0    0    1    1    1    1    1    1    1    1    1    1
##   730558     0    0    0    0    0    0    0    0    1    0    1    1
##   731287     0    0    0    0    0    0    0    0    0    1    1    0
##   732106     0    0    1    0    1    0    0    1    1    1    1    1
##   732107     0    1    0    0    1    0    0    1    1    1    1    1
##   732831     0    0    1    1    1    1    1    0    1    1    1    1
##   732835     1    0    0    0    1    0    1    1    0    1    1    0
##   732836     0    0    0    0    0    1    1    0    1    1    1    1
##   732837     0    0    0    0    0    1    1    1    0    0    0    1
##   733560     0    0    0    1    1    1    0    0    0    0    0    0
##   733565     0    0    0    0    0    0    1    1    1    1    1    1
##   733566     0    0    1    1    0    0    1    1    1    1    1    1
##   733567     0    0    0    0    1    0    0    1    0    1    1    1
##   734211     0    0    0    1    1    1    1    1    1    1    1    1
##   734928     0    0    0    0    0    0    0    0    0    0    1    0
##   734931     0    0    0    0    0    1    0    0    0    1    1    1
##   734940     0    0    0    0    0    0    0    1    1    0    0    0
##   735590     1    0    0    0    0    0    0    0    0    0    0    0
##   735670     1    0    0    0    0    0    0    1    0    0    0    1
##   736387     0    0    0    0    0    0    0    0    0    0    1    0
##   736390     0    0    0    0    0    0    0    0    0    0    1    0
##   736399     0    0    0    0    0    0    1    1    1    1    1    1
##   736400     0    0    0    0    0    0    0    1    1    0    1    1
##   737129     1    0    0    1    1    1    1    1    1    1    1    1
##   737130     0    0    0    0    1    0    1    1    1    1    1    0
##   737782     0    0    0    0    0    0    0    0    1    0    1    0
##   737829     0    0    1    0    0    0    0    0    0    0    0    0
##   737859     0    0    0    1    1    1    1    1    1    1    1    1
##   737860     0    0    0    0    0    0    0    1    1    1    1    0
##   737864     0    1    1    0    1    1    1    1    0    1    1    1
##   738558     0    0    0    0    0    1    1    0    0    0    0    0
##   738589     0    0    0    0    0    0    0    1    0    1    1    0
##   738593     0    0    0    0    0    0    0    0    1    1    1    1
##   738594     0    1    1    1    0    1    1    0    1    1    1    1
##   738668     0    0    1    1    1    1    1    1    1    1    1    1
##   738669     0    0    0    0    1    1    0    1    1    1    1    1
##   738671     0    0    0    0    1    0    0    0    0    1    1    0
##   739289     0    0    0    0    0    0    0    0    1    1    1    1
##   739323     0    0    1    1    1    1    1    1    1    1    1    1
##   739324     0    0    0    0    0    0    0    0    1    1    1    1
##   739397     0    0    0    0    0    0    0    0    1    1    0    1
##   739398     0    0    1    1    0    1    1    0    1    0    0    0
##   739399     0    0    0    0    1    0    0    0    1    1    1    0
##   739400     0    0    0    1    1    1    0    0    0    1    0    0
##   739401     0    0    0    1    1    1    1    0    1    1    1    1
##   739402     1    0    0    0    0    0    0    0    0    0    1    1
##   740017     0    0    0    0    0    0    0    0    1    1    1    0
##   740052     0    0    0    0    1    0    1    1    0    0    1    1
##   740053     0    0    0    0    1    1    1    1    1    1    1    1
##   740127     0    1    0    0    1    1    1    1    1    1    0    1
##   740128     0    1    1    1    1    1    1    1    1    1    0    1
##   740129     0    0    1    0    1    0    0    0    1    0    0    1
##   740130     0    1    1    1    1    1    1    1    1    1    1    1
##   740131     0    0    0    0    1    0    0    0    0    1    0    0
##   740767     0    0    0    0    0    0    0    1    1    1    1    1
##   740768     0    0    0    0    0    0    0    0    1    1    1    1
##   740856     0    1    1    1    0    1    1    1    1    1    1    1
##   740857     0    0    0    0    0    0    0    0    0    1    1    1
##   740858     0    0    0    0    1    0    1    1    1    0    1    1
##   740859     0    1    1    1    1    0    1    1    1    1    1    1
##   740860     0    0    1    1    0    1    0    1    0    1    1    0
##   740861     0    0    0    1    0    0    0    1    0    0    1    0
##   741586     0    0    1    1    0    0    1    1    1    1    1    1
##   741587     0    0    0    0    0    0    0    1    1    0    1    0
##   741588     0    0    0    0    0    1    0    0    1    0    1    1
##   741590     0    0    0    1    1    1    0    0    0    0    1    1
##   742214     0    0    0    0    0    0    0    0    1    0    1    0
##   742227     0    0    0    0    0    0    0    1    1    1    1    0
##   742315     0    0    0    1    0    0    1    0    0    0    0    0
##   742320     0    0    1    1    1    0    1    1    1    0    1    1
##   742321     0    0    0    0    1    1    0    1    0    0    1    0
##   742322     0    0    1    1    1    1    1    0    1    1    0    0
##   743044     0    0    0    0    0    0    1    1    1    1    0    0
##   743045     0    0    0    0    1    0    1    0    0    0    0    1
##   743050     0    0    1    1    1    1    1    1    1    1    1    0
##   743051     0    1    1    1    0    0    0    0    1    1    1    0
##   743052     0    0    0    1    1    1    1    1    1    1    0    0
##   743777     0    0    0    0    0    0    1    1    1    1    1    1
##   743780     0    1    1    1    1    1    1    1    1    1    1    1
##   744506     0    0    0    0    0    0    1    1    1    0    1    1
##   744507     0    0    0    0    0    0    1    1    1    1    1    0
##   745236     0    0    0    0    0    0    1    1    0    0    0    1
##   745962     0    0    0    1    0    1    0    1    1    0    0    0
##   745963     0    0    0    0    0    0    1    0    1    0    0    0
##   746691     1    0    0    1    1    0    0    0    0    0    0    1
##   746692     0    0    0    1    0    0    1    1    0    0    0    0
##   746693     1    0    0    1    1    1    0    0    0    0    0    0
##   747421     0    0    0    1    0    0    0    0    0    1    1    1
##   747423     1    0    0    1    0    0    0    0    0    1    0    1
##   748150     1    0    0    1    0    0    0    0    1    0    1    1
##   748880     0    0    0    0    1    0    0    0    1    0    1    1
##   748881     0    0    0    1    0    0    0    0    0    0    0    0
##   750260     0    0    0    0    0    0    0    0    0    1    0    0
##   753164     0    0    0    0    0    0    0    0    1    1    0    0
##   753165     0    1    0    0    0    0    0    1    1    1    1    1
##   753900     0    0    0    0    1    0    1    1    1    1    1    1
##   753985     0    0    1    0    0    0    1    0    0    0    1    1
##   753986     0    0    0    0    0    0    0    0    1    1    0    0
##   754607     0    0    0    1    0    1    1    0    1    0    1    1
##   754608     0    0    0    0    1    1    1    1    1    1    1    1
##   754715     0    1    1    0    0    0    1    0    0    0    1    1
##   755444     0    0    0    1    1    0    0    1    1    0    0    0
##   755445     0    0    0    0    0    0    0    1    1    0    1    1
##   755446     0    0    0    0    1    0    0    0    0    0    1    1
##   756173     0    0    0    0    0    1    1    1    1    0    1    1
##   756174     0    0    1    0    1    0    0    1    1    0    1    0
##   756175     0    0    0    0    0    1    0    0    0    0    0    0
##   756903     0    0    0    0    0    0    0    0    0    0    0    1
##   757553     1    1    1    0    1    1    1    1    1    1    1    1
##   757563     0    0    0    0    0    0    1    0    0    0    0    0
##   757634     0    0    0    0    0    0    1    0    0    0    1    1
##   758286     0    0    0    0    0    1    1    1    1    1    1    1
##   758363     0    0    0    0    1    0    1    0    1    0    0    1
##   759023     0    1    0    1    0    1    1    1    1    1    1    1
##   759750     1    0    0    0    0    0    1    1    1    1    1    1
##   759751     1    1    1    0    1    1    1    1    1    1    1    1
##   759753     0    0    0    0    1    1    1    1    1    1    1    1
##   760450     1    0    0    1    0    1    0    0    1    1    1    1
##   760548     0    0    1    1    1    1    1    1    1    0    1    1
##   760549     0    0    0    0    1    1    1    0    1    1    1    1
##   761276     0    0    0    0    0    0    0    1    0    0    1    0
##   761277     0    0    1    0    0    0    0    0    0    1    1    1
##   761278     0    0    0    0    0    1    0    1    0    0    1    1
##   761279     0    0    1    1    0    0    1    1    0    0    1    0
##   762006     1    0    1    1    0    1    1    1    1    1    1    1
##   762007     0    0    1    0    0    0    0    1    1    0    1    1
##   762008     0    0    0    0    1    1    1    1    1    1    1    0
##   762736     0    0    0    1    0    0    0    0    1    1    0    0
##   763424     0    0    0    0    0    0    0    1    1    1    1    0
##   763425     0    0    0    0    1    0    0    1    1    1    0    0
##   764886     0    0    0    0    0    0    0    1    0    0    1    0
##   764919     0    0    0    0    0    0    0    1    0    1    0    1
##   764921     0    0    0    0    0    0    0    1    0    0    0    0
##   765617     0    0    0    0    0    1    1    1    1    1    1    1
##   765649     0    1    0    0    0    0    0    1    1    0    1    1
##   765650     0    0    0    0    0    0    0    0    0    0    1    1
##   766345     0    0    0    0    0    0    0    1    1    1    0    0
##   766346     0    0    0    0    0    0    1    0    0    0    0    1
##   766378     0    0    0    0    1    1    1    1    0    1    1    1
##   766379     0    0    0    0    0    0    1    1    0    0    0    1
##   767073     0    0    0    0    0    0    0    0    0    1    1    0
##   767107     0    0    1    0    0    1    0    1    0    0    0    1
##   768562     1    0    0    0    0    0    1    1    0    0    1    1
##   768564     0    1    1    0    0    0    0    0    0    0    1    0
##   769291     0    0    0    0    0    0    0    0    1    1    1    0
##   769292     1    0    0    1    1    0    1    1    0    0    1    1
##   770017     0    0    0    0    0    0    0    0    0    0    1    1
##   770020     0    0    0    0    1    0    1    0    1    1    0    1
##   770746     0    0    0    0    0    0    1    1    1    1    1    0
##   770747     0    0    0    0    0    0    0    1    1    1    1    1
##   770748     0    0    0    0    0    0    0    0    1    1    1    1
##   771412     0    0    0    0    0    0    0    0    0    0    0    1
##   771474     0    0    0    0    0    0    1    1    1    1    1    1
##   771475     0    1    1    0    0    1    1    1    1    1    1    1
##   771476     0    0    0    0    0    0    0    1    1    1    1    1
##   771477     0    0    0    0    0    0    0    1    0    0    0    0
##   772205     0    1    1    0    0    0    0    1    1    0    1    1
##   772958     0    0    1    1    1    1    0    1    0    1    1    0
##   773687     0    1    0    1    1    1    1    1    1    0    1    1
##   773688     0    0    0    0    0    1    1    1    1    1    0    1
##   774416     0    0    0    0    1    1    0    1    1    0    1    1
##   774417     1    0    0    1    1    1    1    1    1    1    1    1
##   775145     0    1    0    1    1    1    1    1    1    1    1    1
##   775146     1    0    1    1    1    1    1    1    1    1    1    1
##   775147     0    0    0    0    1    0    1    1    1    1    0    1
##   775875     0    0    0    1    0    1    0    0    1    0    1    1
##   775876     0    0    0    0    1    1    1    1    1    1    1    1
##   779475     0    0    0    0    0    0    0    0    0    1    0    0
##   782412     1    1    1    1    1    1    1    1    1    1    1    1
##   782413     1    1    1    1    0    1    1    1    1    1    1    1
##   783150     0    0    1    1    1    1    1    1    1    1    1    1
##   783151     1    1    0    0    1    1    1    1    1    0    1    1
##   783152     0    1    1    1    1    1    1    1    1    1    1    1
##   783872     1    1    1    1    1    1    1    1    1    1    1    1
##   783873     1    1    1    1    1    1    1    1    1    1    1    1
##   783880     0    0    1    0    0    0    0    0    0    0    0    1
##   783881     0    0    0    1    0    0    1    1    0    0    0    0
##   783882     1    1    1    1    1    1    1    1    1    1    1    1
##   783883     0    1    1    1    1    1    1    1    1    1    1    1
##   784603     1    1    1    1    1    1    1    1    1    1    1    1
##   785332     0    1    1    1    1    1    1    1    1    1    1    1
##   789709     0    0    0    1    0    0    0    0    0    0    0    0
##   797711     0    0    0    0    0    0    0    0    0    1    1    1
##Filter cells with more than 6 years of sampling
CellsMore6yrs <- CoeflaSimpAnual %>% 
  group_by(Site.k, Year.t) %>%
  summarise(Lists = n()) %>% 
  group_by(Site.k) %>% 
  summarise(List_Year = n()) %>%
  filter(List_Year >= 6)
## `summarise()` has grouped output by 'Site.k'. You can override using the
## `.groups` argument.
#Select only cells with more than 6 years of data (from 3736 to 2695)
CoeflaSimpAnual = CoeflaSimpAnual[CoeflaSimpAnual$Site.k %in% CellsMore6yrs$Site.k, ]

#Chronological order
CoeflaSimpAnual <- CoeflaSimpAnual[order(CoeflaSimpAnual$Site.k, CoeflaSimpAnual$Year.t),]



ggplot(CoeflaSimpAnual, aes(x = Year.t, group = Site.k))+
  geom_line(aes(y = log(Observed.year.y), color = SpRichness), linetype = "dashed")+
  scale_color_gradient2(low = alpha("#4575b4", 0.95), 
                        mid = alpha("#ffffbf", 0.95), midpoint = 73,
                        high = alpha("#d73027", 0.95))+
  geom_point(aes(y = log(Observed.year.y)), size = 0.1, shape = 2)+
  scale_y_continuous(expand = c(0,0))+
  scale_x_continuous(expand = c(0,0), breaks = c(2012, 2014, 2016, 2018, 2020, 2022))+
  labs(x = "Observation date",
       y = "Log(mean observed)",
       color = "Species richness")+
  theme_classic()+
  theme(legend.position = "bottom")

onlye some sites to represent the process and run the model

PopID<- unique(CoeflaSimpAnual$Site.k) 
#308 cells for conducting the example

#empty list per PopID
Coefla_Parm <- data.frame(cell = character(), 
                          mu.hat = numeric(),
                          theta.hat = numeric(),
                          beta.hat = numeric(),
                          tau.hat = numeric())

for(PopID in unique(CoeflaSimpAnual$Site.k)){
  
  dt_tmp <- CoeflaSimpAnual[CoeflaSimpAnual$Site.k == PopID, ]
  
  Observed.t <- dt_tmp$Observed.year.y
  Time.t <- dt_tmp$Year.t
  #-------- Log-transform the observations to carry all the calculations in this program--------------------#
  log.obs    <- log(Observed.t); 
  #---------------------------------------------------------------------------------------------------------#
  
  #----------------------------------------------------------------------
  #        PARAMETER ESTIMATION, PARAMETRIC BOOTSTRAP AND PREDICTIONS
  #----------------------------------------------------------------------
  # Before doing the calculations, the user has to specify ONLY the following 4 options:
  
  # 1. Do you want to compute the ML estimates or the REML estimates?
  
  method <- "REML" # alternatively, set method <- "ML"
  
  # 2. Do you want to plot the predictions?
  pred.plot <- "FALSE" # Set it to "FALSE" if you do not want to plot the predictions
  
  # 3. Do you want to plot the parametric bootstrap distribution of the estimates?
  pboot.plot <- "FALSE" # Set it to "FALSE" if you do not want to plot the bootstrap distribution of the estimates
  
  # 4. How many bootstrap replicates?
  NBoot <- 10; # 10 just to try, for formal results use 1000 or more 
  
  #-------------------------------------------------------------------------------------------------#
  #  5. RUN THE FOLLOWING LINE OF CODE TO COMPUTE THE ESTIMATES, PREDICTIONS,
  #     AND RUN A PARAMETRIC BOOTSTRAP. THE USER DOES NOT NEED TO MODIFY THIS LINE OF CODE.
  #     THE OUTPUT OF THE FUNCTION 'ROUSS.CALCS' IS A LIST AND THE USER CAN RETRIEVE EACH OF THE   
  #     LIST ELEMENTS PRINTED AND SAVED IN THE OBJECT "all.results". 
  #     THE 95% PARAMETRIC BOOTSTRAP FOR BOTH, THE PARAMETERS AND THE PREDICTIONS ARE COMPUTED 
  #     BY THE FUNCTION "ROUSS.CALCS"   
  #      
  #-------------------------------------------------------------------------------------------------#
  
  all.results <- ROUSS.CALCS(Yobs=log.obs,Tvec=Time.t, pmethod=method, nboot=NBoot, plot.pred=pred.plot, plot.bootdists = pboot.plot)
  
  Coefla_Parm <- rbind(Coefla_Parm, c(cell = PopID, all.results$parms.est))
  
  print(all.results$parms.est)
}

#It changes the column name, but save the estimates!!
colnames(Coefla_Parm) = c("cell", "mu.hat", "theta.hat", "beta.hat", "tau.hat")

saveRDS(Coefla_Parm, "Coefla_ParametersOUSS_abundance.rds")

\(\hat{\mu}\) is the mean stationary log-abundance estimate for each site!

We can also predict the estimated abundance for each year:

PopID<- unique(CoeflaSimpAnual$Site.k) 

GSS.Ban <- list()

for(PopID in unique(CoeflaSimpAnual$Site.k)){
  
  ParmPop <- Coefla_Parm %>%
    filter(cell == PopID) %>%
    dplyr::select(mu.hat, theta.hat, beta.hat, tau.hat)

    dt_tmp <- CoeflaSimpAnual[CoeflaSimpAnual$Site.k == PopID, ]
  
    Observed.t <- dt_tmp$Observed.year.y
    Tvec <- dt_tmp$Year.t
  #-------- Log-transform the observations to carry all the calculations in this program--------------------#
    Yobs <- log(Observed.t); 

    ROUSS.predict(parms = as.numeric(ParmPop), Yobs = Yobs, Tvec = Tvec, plot.it = "TRUE")

  PRBan.parms <- c(as.numeric(ParmPop))
  #Parameteric Bootstrap function and confidence intervals
  B = NBoot

  tt <- Tvec-Tvec[1]
  long.t <- tt[1]:max(tt)
  nparms    <- length(PRBan.parms);
  preds.boot1<- matrix(0,nrow=B,ncol=length(tt))
  preds.boot2<- matrix(0,nrow=B,ncol=length(long.t))

  boot.remles <- matrix(0,nrow=B,ncol=nparms+1); 
  all.sims  <- ROUSS.sim(nsims=B,parms=PRBan.parms,Tvec=Tvec);
  all.preds <- ROUSS.predict(parms=PRBan.parms, Yobs=Yobs,Tvec=Tvec,plot.it="FALSE")
  reml.preds<- all.preds[[1]][,2]
  reml.longpreds <- all.preds[[2]][,2]      

  for(b in 1:B ){
  
    bth.timeseries <- all.sims[,b];
    remles.out <- ROUSS.REML(Yobs=bth.timeseries,Tvec=Tvec, parms.guess=PRBan.parms);
    boot.remles[b,] <- c(remles.out$remls, remles.out$lnLhat);
    all.bootpreds <- ROUSS.predict(parms=remles.out$remls, Yobs=bth.timeseries,Tvec=Tvec,plot.it="FALSE");
    preds.boot1[b,] <- all.bootpreds[[1]][,2]
    preds.boot2[b,] <- all.bootpreds[[2]][,2]
  
  } 

CIs.mat <- apply(boot.remles,2,FUN=function(x){quantile(x,probs=c(0.025,0.975))});
CIs.mat <- rbind(CIs.mat[1,1:4],PRBan.parms,CIs.mat[2,1:4]);
rownames(CIs.mat) <- c("2.5%","REMLE","97.5%");
colnames(CIs.mat) <- c("mu", "theta","betasq","tausq");

preds.CIs1 <- apply(preds.boot1,2,FUN=function(x){quantile(x,probs=c(0.025,0.975))});
mean.boots <- apply(preds.boot1,2,FUN=function(x){quantile(x,probs=0.50)})
preds.CIs1 <- t(rbind(Tvec,reml.preds-(mean.boots-preds.CIs1[1,]), reml.preds, reml.preds+(preds.CIs1[2,]-mean.boots)));
colnames(preds.CIs1) <- c("Year","2.5%","REMLE","97.5%");

preds.CIs2 <- apply(preds.boot2,2,FUN=function(x){quantile(x,probs=c(0.025,0.975))});
mean.boots2 <- apply(preds.boot2,2,FUN=function(x){quantile(x,probs=0.50)})
preds.CIs2 <- t(rbind(long.t+Tvec[1],reml.longpreds-(mean.boots2-preds.CIs2[1,]), reml.longpreds, reml.longpreds+(preds.CIs2[2,]-mean.boots2)));

#preds.CIs2 <- apply(preds.boot2,2,FUN=function(x){quantile(x,probs=c(0.025,0.5,0.975))});
#preds.CIs2 <- t(rbind(long.t+Tvec[1],preds.CIs2));
colnames(preds.CIs2) <- c("Year","2.5%","REMLE","97.5%");
#pred.CIs2 <- cbind(preds.CIs2[,1], reml.longpreds-(preds.CIs2[,3]-preds.CIs2[,2]),reml.longpreds,reml.longpreds+(preds.CIs2[,4]-preds.CIs2[,3]))

boot.list <- list(boot.remles = boot.remles, CIs.mat = CIs.mat, preds.CIs1 = preds.CIs1, 
                  preds.CIs2=preds.CIs2)


boot.list$preds.CIs1

df.boot <- as.data.frame(cbind(PopID, boot.list$preds.CIs1))

GSS.Ban <- c(GSS.Ban, list(df.boot))

}


#extract estimates
GSS.BanAll <- bind_rows(GSS.Ban)
colnames(GSS.BanAll) <-c("Site.k", "Year.t","Low","REMLE","High");

GSS.BanAll <- GSS.BanAll %>%
  mutate_at(c("Year.t","Low","REMLE","High"), as.numeric) %>%
  left_join(CoeflaSimpAnual, by = c("Site.k", "Year.t"))

saveRDS(GSS.BanAll, "Coefla_predictedAbundance_REMLE_OUSS.rds")
GSS.BanAll <- readRDS("Coefla_predictedAbundance_REMLE_OUSS.rds")

ggplot(GSS.BanAll, aes(x = Year.t))+
  geom_line(aes(y = log(Observed.year.y), group = Site.k, color = SpRichness), linetype = "dotted")+
  geom_point(aes(y = log(Observed.year.y), fill = SpRichness), shape = 21)+
  geom_line(aes(y = log(REMLE), group = Site.k, color = SpRichness), linetype = "solid")+
  geom_point(aes(y = log(REMLE), fill = SpRichness), shape = 22)+
  scale_color_gradient2(low = alpha("#4575b4", 0.95), 
                        mid = alpha("#ffffbf", 0.95), midpoint = 73,
                        high = alpha("#d73027", 0.95))+
  scale_fill_gradient2(low = alpha("#4575b4", 0.95), 
                        mid = alpha("#ffffbf", 0.95), midpoint = 73,
                        high = alpha("#d73027", 0.95))+
  labs(x = "t",
       y = "log(Observed / estimated abundance)",
       color = "Species richness",
       fill = "Species richness",
       title = "62 sites (hexagonal cells)")+
  theme_classic()+
  theme(legend.position = c(0.8,0.8),
        legend.direction = "horizontal")
## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

ggplot(GSS.BanAll, aes(x = SpRichness, fill = SpRichness))+
  geom_point(aes(y = REMLE), shape = 22)+
  geom_point(aes(y = Observed.year.y), shape = 21)+
  geom_smooth(aes(x = SpRichness, y = REMLE), method = "lm", linetype = "solid")+
  geom_smooth(aes(x = SpRichness, y = Observed.year.y), method = "lm", linetype = "dashed")+
  scale_fill_gradient2(low = alpha("#4575b4", 0.95), 
                        mid = alpha("#ffffbf", 0.95), midpoint = 73,
                        high = alpha("#d73027", 0.95))+
  labs(x = "Species richness",
       y = "Observed / estimated abundance",
       fill = "Species richness",
       title = "Relationship of abundance and species richness")+
  theme_classic()+
  theme(legend.position = c(0.75,0.75),
        legend.direction = "horizontal")
## `geom_smooth()` using formula = 'y ~ x'
## Warning: The following aesthetics were dropped during statistical transformation: fill.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?
## `geom_smooth()` using formula = 'y ~ x'
## Warning: The following aesthetics were dropped during statistical transformation: fill.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

We can map our data

GSS.map.df <- GSS.BanAll %>%
  group_by(Site.k, Year.t, REMLE) %>%
  summarise(SpRichness = mean(SpRichness),
            lon = mean(longitude),
            lat = mean(latitude))
## `summarise()` has grouped output by 'Site.k', 'Year.t'. You can override using
## the `.groups` argument.
ggplot() +
  geom_sf(data = world1)+
  geom_point(data = GSS.map.df, 
             aes(x = lon, y = lat, size = REMLE, fill = SpRichness), alpha = 0.5, shape = 21)+
  scale_fill_gradient2(low = alpha("#4575b4", 0.75), 
                       mid = alpha("#ffffbf", 0.75), midpoint = 73,
                       high = alpha("#d73027", 0.75))+
  coord_sf(xlim = c(-91, -59.4), 
           ylim =  c(8, 27)) +
  scale_x_continuous(breaks = seq(from = -90, to = -60, by = 15))+
  labs(y = "Latitude",
       x = "Longitude",       
       tag = "",
       title = "Caribbean (high eBird coverage-based completeness)",
       subtitle = "Abundance estimation using GSS(OUSS) models for Bananaquit",
       fill = "Species richness",
       size = "Stationary abundance estimation") +
  theme_classic()+
  theme(legend.position = "bottom",
        legend.direction = "horizontal")+
  geom_magnify(data = GSS.map.df, aes(x = lon, y = lat), 
               from = c(xmin = -86.2, xmax = -81.2, ymin = 8, ymax = 13.5), 
               to = c(xmin = -75, xmax = -62, ymin = 10, ymax = 25),
               shadow = TRUE) 
## Warning in layer_sf(geom = ggproto(NULL, GeomMagnify), mapping = mapping, :
## Ignoring unknown aesthetics: x and y
## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

## Warning in .setMask(mask$f, mask$ref): Ignored alpha mask (not supported on
## this device)

End of this document. Orlando.

A previous approach (reducing the observer error parameter and applying Data Cloning)

We tried with the Poisson process (adjust to the length of the time series). Here is the GSS.dc model:

StochGSS2.dc <- function(){

  # Priors on model parameters. Priors are DC1 in Lele et al (2007)
    a1 ~ dnorm(0,1);   # constant, the population growth rate. 
    c1 ~ dunif(-1,1);      # constant, the density dependence parameter. 
    sig1 ~ dlnorm(-0.5,10); #variance parameter of stochastic  environment (process noise) in the system
    stovar1 <- 1/pow(sig1,2)

    for(k in 1:K){
      # Simulate trajectory that depends on the previous
      mean_X1[1,k] <- a1/(1-c1) # Expected value of the first realization of the process
                          # this is drawn from the stationary distribution of the process
                          # Equation 14 (main text) and  A.4 in Appendix of Dennis et al 2006
      Varno1[k] <- pow(sig1,2)/(1-pow(c1,2)) #. Equation A.5 in Appendix of Dennis et al 2006
      
 
      # Updating the state: Stochastic process for 29 steps
      X1[1,k]~dnorm(mean_X1[1,k], 1/Varno1[k]); #first estimation of population
      
      #iteration of the GSS model in the data
      for (i in 2:12) {
        mean_X1[i,k] <- a1 + c1 * X1[(i - 1),k]
        X1[i,k] ~ dnorm(mean_X1[i,k], stovar1) # Et is included here since a+cX(t-1) + Et ~ Normal(a+cX(t-1),sigma^2)
      }
  
    # Updating the observations, from the counts
      for (i in 1:12) {
        Y1[i,k] ~ dpois(exp(X1[i,k])) #Here we can plug the probability of detection, p*exp(X1[i,k]), and we could also estimate it from the data!
      }
  }

}

Now, we select the data, first Costa Rica - Alejuela (CR-A), which has complete data 2012-2023.

table(Pop$state_code, Pop$year)

#Costa Rica - Alejuela

CR.A <- Pop %>%
  dplyr::filter(state_code == "CR-A") %>%
  dplyr::select(mu_count)

CR <- round(CR.A$mu_count, 0)

datalistGSS.CR.dc <- list(K = 1,
                       Y1 = dcdim(data.matrix(CR))) # Not need to be in LOG!!!

dcrun.GSS2.CR <- dc.fit(data = datalistGSS.CR.dc,
                          params = c("a1", "c1", "sig1"), #extract from the model
                          model = StochGSS2.dc, 
                          n.clones = c(5,10,50,100,150),
                          multiply = "K",
                          n.chains = 10,
                          n.adapt = 10000,
                          n.update = 10000,
                          thin = 10,
                          n.iter = 100000)

summary(dcrun.GSS2.CR);

dcdiag(dcrun.GSS2.CR) 

Now, Aruba (Lesser Antilles), missing the two first years of data, complete data 2014-2023

A <- Pop %>%
  dplyr::filter(state_code == "AW-") %>%
  dplyr::select(mu_count)

AW <- c(round(A$mu_count, 0)) #the

#and have to adjust the model
StochGSS2.AW.dc <- function(){

  # Priors on model parameters. Priors are DC1 in Lele et al (2007)
    a1 ~ dnorm(0,1);   # constant, the population growth rate. 
    c1 ~ dunif(-1,1);      # constant, the density dependence parameter. 
    sig1 ~ dlnorm(-0.5,10); #variance parameter of stochastic  environment (process noise) in the system
    stovar1 <- 1/pow(sig1,2)

    for(k in 1:K){
      # Simulate trajectory that depends on the previous
      mean_X1[1,k] <- a1/(1-c1) # Expected value of the first realization of the process
                          # this is drawn from the stationary distribution of the process
                          # Equation 14 (main text) and  A.4 in Appendix of Dennis et al 2006
      Varno1[k] <- pow(sig1,2)/(1-pow(c1,2)) #. Equation A.5 in Appendix of Dennis et al 2006
      
 
      # Updating the state: Stochastic process for 29 steps
      X1[1,k]~dnorm(mean_X1[1,k], 1/Varno1[k]); #first estimation of population
      
      #iteration of the GSS model in the data
      for (i in 2:10) {
        mean_X1[i,k] <- a1 + c1 * X1[(i - 1),k]
        X1[i,k] ~ dnorm(mean_X1[i,k], stovar1) # Et is included here since a+cX(t-1) + Et ~ Normal(a+cX(t-1),sigma^2)
      }
  
    # Updating the observations, from the counts
      for (i in 1:10) {
        Y1[i,k] ~ dpois(exp(X1[i,k])) #Here we can plug the probability of detection, p*exp(X1[i,k]), and we could also estimate it from the data!
      }
  }

}


datalistGSS.AW.dc <- list(K = 1,
                       Y1 = dcdim(data.matrix(AW))) # Not need to be in LOG!!!

dcrun.GSS2.AW <- dc.fit(data = datalistGSS.AW.dc,
                          params = c("a1", "c1", "sig1"), #extract from the model
                          model = StochGSS2.AW.dc, 
                          n.clones = c(5,10,50,100,150),
                          multiply = "K",
                          n.chains = 10,
                          n.adapt = 10000,
                          n.update = 10000,
                          thin = 10,
                          n.iter = 100000)

summary(dcrun.GSS2.AW);

dcdiag(dcrun.GSS2.AW) 

Let’s try Puerto Rico - San Juan, which has NA in 2013 and 2014. I will run from 2015 to 2023

PR.SJ <- Pop %>%
  dplyr::filter(state_code == "PR-SJ") %>%
  dplyr::select(mu_count)

PR <- c(round(PR.SJ$mu_count, 0)) #the

#and have to adjust the model
StochGSS2.PR.dc <- function(){

  # Priors on model parameters. Priors are DC1 in Lele et al (2007)
    a1 ~ dnorm(0,1);   # constant, the population growth rate. 
    c1 ~ dunif(-1,1);      # constant, the density dependence parameter. 
    sig1 ~ dlnorm(-0.5,10); #variance parameter of stochastic  environment (process noise) in the system
    stovar1 <- 1/pow(sig1,2)

    for(k in 1:K){
      # Simulate trajectory that depends on the previous
      mean_X1[1,k] <- a1/(1-c1) # Expected value of the first realization of the process
                          # this is drawn from the stationary distribution of the process
                          # Equation 14 (main text) and  A.4 in Appendix of Dennis et al 2006
      Varno1[k] <- pow(sig1,2)/(1-pow(c1,2)) #. Equation A.5 in Appendix of Dennis et al 2006
      
 
      # Updating the state: Stochastic process for 29 steps
      X1[1,k]~dnorm(mean_X1[1,k], 1/Varno1[k]); #first estimation of population
      
      #iteration of the GSS model in the data
      for (i in 2:10) {
        mean_X1[i,k] <- a1 + c1 * X1[(i - 1),k]
        X1[i,k] ~ dnorm(mean_X1[i,k], stovar1) # Et is included here since a+cX(t-1) + Et ~ Normal(a+cX(t-1),sigma^2)
      }
  
    # Updating the observations, from the counts
      for (i in 1:10) {
        Y1[i,k] ~ dpois(exp(X1[i,k])) #Here we can plug the probability of detection, p*exp(X1[i,k]), and we could also estimate it from the data!
      }
  }

}


datalistGSS.PR.dc <- list(K = 1,
                       Y1 = dcdim(data.matrix(PR))) # Not need to be in LOG!!!

dcrun.GSS2.PR <- dc.fit(data = datalistGSS.PR.dc,
                          params = c("a1", "c1", "sig1"), #extract from the model
                          model = StochGSS2.PR.dc, 
                          n.clones = c(5,10,50,100,150),
                          multiply = "K",
                          n.chains = 10,
                          n.adapt = 10000,
                          n.update = 10000,
                          thin = 10,
                          n.iter = 100000)

summary(dcrun.GSS2.PR);

dcdiag(dcrun.GSS2.PR) 

Let’s try Venezuela - Miranda, which has data from 2016 to 2023

VE.M <- Pop %>%
  dplyr::filter(state_code == "VE-M") %>%
  dplyr::select(mu_count)

VE <- c(round(VE.M$mu_count, 0)) #the

#and have to adjust the model
StochGSS2.VE.dc <- function(){

  # Priors on model parameters. Priors are DC1 in Lele et al (2007)
    a1 ~ dnorm(0,1);   # constant, the population growth rate. 
    c1 ~ dunif(-1,1);      # constant, the density dependence parameter. 
    sig1 ~ dlnorm(-0.5,10); #variance parameter of stochastic  environment (process noise) in the system
    stovar1 <- 1/pow(sig1,2)

    for(k in 1:K){
      # Simulate trajectory that depends on the previous
      mean_X1[1,k] <- a1/(1-c1) # Expected value of the first realization of the process
                          # this is drawn from the stationary distribution of the process
                          # Equation 14 (main text) and  A.4 in Appendix of Dennis et al 2006
      Varno1[k] <- pow(sig1,2)/(1-pow(c1,2)) #. Equation A.5 in Appendix of Dennis et al 2006
      
 
      # Updating the state: Stochastic process for 29 steps
      X1[1,k]~dnorm(mean_X1[1,k], 1/Varno1[k]); #first estimation of population
      
      #iteration of the GSS model in the data
      for (i in 2:9) {
        mean_X1[i,k] <- a1 + c1 * X1[(i - 1),k]
        X1[i,k] ~ dnorm(mean_X1[i,k], stovar1) # Et is included here since a+cX(t-1) + Et ~ Normal(a+cX(t-1),sigma^2)
      }
  
    # Updating the observations, from the counts
      for (i in 1:9) {
        Y1[i,k] ~ dpois(exp(X1[i,k])) #Here we can plug the probability of detection, p*exp(X1[i,k]), and we could also estimate it from the data!
      }
  }

}


datalistGSS.VE.dc <- list(K = 1,
                       Y1 = dcdim(data.matrix(VE))) # Not need to be in LOG!!!

dcrun.GSS2.VE <- dc.fit(data = datalistGSS.VE.dc,
                          params = c("a1", "c1", "sig1"), #extract from the model
                          model = StochGSS2.VE.dc, 
                          n.clones = c(5,10,50,100,150),
                          multiply = "K",
                          n.chains = 10,
                          n.adapt = 10000,
                          n.update = 10000,
                          thin = 10,
                          n.iter = 100000)

summary(dcrun.GSS2.VE);
dcdiag(dcrun.GSS2.VE)

Let’s try Grenada - Saint Andrew, which has six years of data

GR.SA <- Pop %>%
  dplyr::filter(state_code == "GD-01") %>%
  dplyr::select(mu_count)

GR <- c(round(GR.SA$mu_count, 0)) #the

#and have to adjust the model
StochGSS2.GR.dc <- function(){

  # Priors on model parameters. Priors are DC1 in Lele et al (2007)
    a1 ~ dnorm(0,1);   # constant, the population growth rate. 
    c1 ~ dunif(-1,1);      # constant, the density dependence parameter. 
    sig1 ~ dlnorm(-0.5,10); #variance parameter of stochastic  environment (process noise) in the system
    stovar1 <- 1/pow(sig1,2)

    for(k in 1:K){
      # Simulate trajectory that depends on the previous
      mean_X1[1,k] <- a1/(1-c1) # Expected value of the first realization of the process
                          # this is drawn from the stationary distribution of the process
                          # Equation 14 (main text) and  A.4 in Appendix of Dennis et al 2006
      Varno1[k] <- pow(sig1,2)/(1-pow(c1,2)) #. Equation A.5 in Appendix of Dennis et al 2006
      
 
      # Updating the state: Stochastic process for 29 steps
      X1[1,k]~dnorm(mean_X1[1,k], 1/Varno1[k]); #first estimation of population
      
      #iteration of the GSS model in the data
      for (i in 2:6) {
        mean_X1[i,k] <- a1 + c1 * X1[(i - 1),k]
        X1[i,k] ~ dnorm(mean_X1[i,k], stovar1) # Et is included here since a+cX(t-1) + Et ~ Normal(a+cX(t-1),sigma^2)
      }
  
    # Updating the observations, from the counts
      for (i in 1:6) {
        Y1[i,k] ~ dpois(exp(X1[i,k])) #Here we can plug the probability of detection, p*exp(X1[i,k]), and we could also estimate it from the data!
      }
  }

}


datalistGSS.GR.dc <- list(K = 1,
                       Y1 = dcdim(data.matrix(GR))) # Not need to be in LOG!!!

dcrun.GSS2.GR <- dc.fit(data = datalistGSS.GR.dc,
                          params = c("a1", "c1", "sig1"), #extract from the model
                          model = StochGSS2.GR.dc, 
                          n.clones = c(5,10,50,100,150),
                          multiply = "K",
                          n.chains = 10,
                          n.adapt = 10000,
                          n.update = 10000,
                          thin = 10,
                          n.iter = 100000)

summary(dcrun.GSS2.GR);

dcdiag(dcrun.GSS2.GR)

And Jamaica - Saint Andrew, which has four years of data

JM.SA <- Pop %>%
  dplyr::filter(state_code == "JM-02") %>%
  dplyr::select(mu_count)

JM <- c(round(JM.SA$mu_count, 0)) #the

#and have to adjust the model
StochGSS2.JM.dc <- function(){

  # Priors on model parameters. Priors are DC1 in Lele et al (2007)
    a1 ~ dnorm(0,1);   # constant, the population growth rate. 
    c1 ~ dunif(-1,1);      # constant, the density dependence parameter. 
    sig1 ~ dlnorm(-0.5,10); #variance parameter of stochastic  environment (process noise) in the system
    stovar1 <- 1/pow(sig1,2)

    for(k in 1:K){
      # Simulate trajectory that depends on the previous
      mean_X1[1,k] <- a1/(1-c1) # Expected value of the first realization of the process
                          # this is drawn from the stationary distribution of the process
                          # Equation 14 (main text) and  A.4 in Appendix of Dennis et al 2006
      Varno1[k] <- pow(sig1,2)/(1-pow(c1,2)) #. Equation A.5 in Appendix of Dennis et al 2006
      
 
      # Updating the state: Stochastic process for 29 steps
      X1[1,k]~dnorm(mean_X1[1,k], 1/Varno1[k]); #first estimation of population
      
      #iteration of the GSS model in the data
      for (i in 2:4) {
        mean_X1[i,k] <- a1 + c1 * X1[(i - 1),k]
        X1[i,k] ~ dnorm(mean_X1[i,k], stovar1) # Et is included here since a+cX(t-1) + Et ~ Normal(a+cX(t-1),sigma^2)
      }
  
    # Updating the observations, from the counts
      for (i in 1:4) {
        Y1[i,k] ~ dpois(exp(X1[i,k])) #Here we can plug the probability of detection, p*exp(X1[i,k]), and we could also estimate it from the data!
      }
  }

}


datalistGSS.JM.dc <- list(K = 1,
                       Y1 = dcdim(data.matrix(JM))) # Not need to be in LOG!!!

dcrun.GSS2.JM <- dc.fit(data = datalistGSS.JM.dc,
                          params = c("a1", "c1", "sig1"), #extract from the model
                          model = StochGSS2.JM.dc, 
                          n.clones = c(5,10,50,100,150),
                          multiply = "K",
                          n.chains = 10,
                          n.adapt = 10000,
                          n.update = 10000,
                          thin = 10,
                          n.iter = 100000)

summary(dcrun.GSS2.JM);

dcdiag(dcrun.GSS2.JM)

Now the null:

H0 <- round(Pop$mu_count, 0)

#and have to adjust the model
StochGSS2.H0.dc <- function(){

  # Priors on model parameters. Priors are DC1 in Lele et al (2007)
    a1 ~ dnorm(0,1);   # constant, the population growth rate. 
    c1 ~ dunif(-1,1);      # constant, the density dependence parameter. 
    sig1 ~ dlnorm(-0.5,10); #variance parameter of stochastic  environment (process noise) in the system
    stovar1 <- 1/pow(sig1,2)

    for(k in 1:K){
      # Simulate trajectory that depends on the previous
      mean_X1[1,k] <- a1/(1-c1) # Expected value of the first realization of the process
                          # this is drawn from the stationary distribution of the process
                          # Equation 14 (main text) and  A.4 in Appendix of Dennis et al 2006
      Varno1[k] <- pow(sig1,2)/(1-pow(c1,2)) #. Equation A.5 in Appendix of Dennis et al 2006
      
 
      # Updating the state: Stochastic process for 29 steps
      X1[1,k]~dnorm(mean_X1[1,k], 1/Varno1[k]); #first estimation of population
      
      #iteration of the GSS model in the data
      for (i in 2:51) {
        mean_X1[i,k] <- a1 + c1 * X1[(i - 1),k]
        X1[i,k] ~ dnorm(mean_X1[i,k], stovar1) # Et is included here since a+cX(t-1) + Et ~ Normal(a+cX(t-1),sigma^2)
      }
  
    # Updating the observations, from the counts
      for (i in 1:51) {
        Y1[i,k] ~ dpois(exp(X1[i,k])) #Here we can plug the probability of detection, p*exp(X1[i,k]), and we could also estimate it from the data!
      }
  }

}


datalistGSS.H0.dc <- list(K = 1,
                       Y1 = dcdim(data.matrix(H0))) # Not need to be in LOG!!!

dcrun.GSS2.H0 <- dc.fit(data = datalistGSS.H0.dc,
                          params = c("a1", "c1", "sig1"), #extract from the model
                          model = StochGSS2.H0.dc, 
                          n.clones = c(5,10,50,100,150),
                          multiply = "K",
                          n.chains = 10,
                          n.adapt = 10000,
                          n.update = 10000,
                          thin = 10,
                          n.iter = 100000)

summary(dcrun.GSS2.H0);

dcdiag(dcrun.GSS2.H0)

Let’s combine the results

MLE <- as.data.frame(rbind(coef(dcrun.GSS2.H0),
             coef(dcrun.GSS2.CR),
             coef(dcrun.GSS2.VE),
             coef(dcrun.GSS2.AW),
             coef(dcrun.GSS2.GR),
             coef(dcrun.GSS2.PR),
             coef(dcrun.GSS2.JM)))

MLE$Pop <- c("Null",
             "CR-A",
             "VE-M",
             "AW-",
             "GD-01",
             "PR-SJ",
             "JM-02")

MLE$N_equ <- exp(-(MLE$a1/(MLE$c1-1)))

print(MLE)

#Some figures

#how many records per site?
Pop %>% 
  group_by(state_code) %>% 
  summarise(AllRecords = sum(n_records))

Mainland sites

par(mfrow = c(1,2));#create a 2x2 plotting matrix

hist(rnorm(10000, mean = MLE[2,5], sd = MLE[2,3]/(1-MLE[2,2])),
     xlab = "Equilibrium density",
     main = "Histogram of Costa Rica - Alejuela",
     xlim = c(0,10));
abline(v=MLE[2,5], col="blue");
abline(v=mean(CR), col = "red");
legend(5, 1500, legend=c("Abundance (estimated)", "Observation mean"),
       col=c("blue", "red"), lty=1, cex=0.8);
text(x = 8, y = 500, "n = 249")

hist(rnorm(10000, mean = MLE[3,5], sd = MLE[3,3]/(1-MLE[3,2])),
     xlab = "Equilibrium density",
     main = "Histogram of Venezuela - Miranda",
     xlim = c(0,10));
abline(v=MLE[3,5], col="blue");
abline(v=mean(VE), col = "red");
text(x = 8, y = 500, "n = 521")

Lesser Antilles sites

par(mfrow = c(1,2));#create a 2x2 plotting matrix

hist(rnorm(10000, mean = MLE[4,5], sd = MLE[4,3]/(1-MLE[4,2])),
     xlab = "Equilibrium density",
     main = "Histogram of Aruba",
     xlim = c(0,10));
abline(v=MLE[4,5], col="blue");
abline(v=mean(AW), col = "red");
text(x = 8, y = 500, "n = 229")

hist(rnorm(10000, mean = MLE[5,5], sd = MLE[5,3]/(1-MLE[5,2])),
     xlab = "Equilibrium density",
     main = "Histogram of Grenada - Saint Andrew",
     xlim = c(0,10));
abline(v=MLE[5,5], col="blue");
abline(v=mean(GR), col = "red");
text(x = 8, y = 500, "n = 117")

Greater Antilles sites

par(mfrow = c(1,2));#create a 2x2 plotting matrix

hist(rnorm(10000, mean = MLE[6,5], sd = MLE[6,3]/(1-MLE[6,2])),
     xlab = "Equilibrium density",
     main = "Histogram of Puerto Rico - San Juan",
     xlim = c(0,10));
abline(v=MLE[6,5], col="blue");
abline(v=mean(PR), col = "red");
text(x = 8, y = 500, "n = 105")

hist(rnorm(10000, mean = MLE[7,5], sd = MLE[7,3]/(1-MLE[7,2])),
     xlab = "Equilibrium density",
     main = "Histogram of Jamaica - Saint Andrew",
     xlim = c(0,10));
abline(v=MLE[7,5], col="blue");
abline(v=mean(JM), col = "red");
text(x = 8, y = 500, "n = 18")

Null hypothesis and relationship

par(mfrow = c(1,2));#create a 2x2 plotting matrix

hist(rnorm(10000, mean = MLE[1,5], sd = MLE[1,3]/(1-MLE[1,2])),
     xlab = "Equilibrium density",
     main = "Histogram of all the Meta-archipelago of the Caribbean",
     xlim = c(0,10));
abline(v=MLE[1,5], col="blue");
abline(v=mean(H0), col = "red");
text(x = 8, y = 500, "n = 1,239")

MLE$S <- c(NA, 149, 76, 67, 51, 39, 31)

plot(x = MLE[c(2:7),6], y = MLE[c(2:7),5], pch = 19, cex = 4,
     xlim = c(0,150), ylim = c(0,5),
     xlab = "Species richness",
     ylab = "Abundance (estimated)",
     main = "Hypothetical relationship",
     col = c("#F36F20", "#DD692E", "#864268", "#704C76", "#434174","#01458A"));
abline(lm(MLE$N_equ ~ MLE$S), col = "darkgray");
text(x = 25, y = 1, expression('R' ^ 2~"= 0.189; p = 0.215"))

Likelihood Ratio test

The parameters for the estimation of the unknown, mean number of Bananaquit per list of ≤2 km and 10-120 minutes in the six sites \(s\) (\(s = CR, VE, AR, GD, PR, JM\)) are \(\mathbf{\Theta} = [a, c, \sigma]\). Several possible scenarios regarding the comparison of these quantities across sites comes to mind. The simple scenario is that all parameters are equal among the meta-archipelago and representing identical parameters \(\mathbf{\Theta}\), thus our working null hypothesis is

\[ {\rm H_{0}} : \mathbf{\Theta}_{CR} = \mathbf{\Theta}_{VE} = \mathbf{\Theta}_{AW} = \mathbf{\Theta}_{GD} = \mathbf{\Theta}_{PR} = \mathbf{\Theta}_{JM} = \mathbf{\Theta} \]

A very contrasting alternative hypothesis (Model1) represent differences among the sites \(s\) in the meta-archipelago, and is \[ {\rm H_{1}} : \mathbf{\Theta}_{CR} \neq \mathbf{\Theta}_{VE} \neq \mathbf{\Theta}_{AW} \neq \mathbf{\Theta}_{GD} \neq \mathbf{\Theta}_{PR} \neq \mathbf{\Theta}_{JM} \]

A third hypothesis is that sub-regions within the meta-archipelago (Mainland vs Lesser Antilles vs Greater Antilles) have similar parameters \(\mathbf{\Theta}_{Subregion}\), and they differ from the others…

\[ {\rm H_{2}} : \mathbf{\Theta}_{CR} = \mathbf{\Theta}_{VE} \neq \mathbf{\Theta}_{AW} = \mathbf{\Theta}_{GD} \neq \mathbf{\Theta}_{PR} = \mathbf{\Theta}_{JM} \] or \[ {\rm H_{2}} : \mathbf{\Theta}_{Mainland} \neq \mathbf{\Theta}_{Lesser Antilles} \neq \mathbf{\Theta}_{Greater Antilles} \]

And another hypothesis is that representing “the island syndrome” (Mainland vs Islands) of difference in parameters \(\mathbf{\Theta}_{Mainland-Island}\):

\[ {\rm H_{3}} : \mathbf{\Theta}_{CR} = \mathbf{\Theta}_{VE} \neq \mathbf{\Theta}_{AW} = \mathbf{\Theta}_{GD} = \mathbf{\Theta}_{PR} = \mathbf{\Theta}_{JM} \] or \[ {\rm H_{3}} : \mathbf{\Theta}_{Mainland} \neq \mathbf{\Theta}_{Islands} \]

The LR test statistics are \[ {\Lambda_{0-1}} = \frac{L_0(\mathbf{\hat\Theta})}{L_1(\mathbf{\hat\Theta}_{CR},\mathbf{\hat\Theta}_{VE},\mathbf{\hat\Theta}_{AW},\mathbf{\hat\Theta}_{GD},\mathbf{\hat\Theta}_{PR},\mathbf{\hat\Theta}_{JM})} \]

\[ {\Lambda_{0-2}} = \frac{L_0(\hat\Theta)}{L_2(\mathbf{\hat\Theta}_{Mainland},\mathbf{\hat\Theta}_{Lesser Antilles},\mathbf{\hat\Theta}_{Greater Antilles})} \]

\[ {\Lambda_{0-3}} = \frac{L_0(\hat\Theta)}{L_3(\mathbf{\hat\Theta}_{Mainland},\mathbf{\hat\Theta}_{Islands})} \]

\[ {\Lambda_{2-1}} = \frac{L_2(\mathbf{\hat\Theta}_{Mainland},\mathbf{\hat\Theta}_{Lesser Antilles},\mathbf{\hat\Theta}_{Greater Antilles})}{L_1(\mathbf{\hat\Theta}_{CR},\mathbf{\hat\Theta}_{VE},\mathbf{\hat\Theta}_{AW},\mathbf{\hat\Theta}_{GD},\mathbf{\hat\Theta}_{PR},\mathbf{\hat\Theta}_{JM})} \]